library(tidyverse)
library(zoo)
library(survminer)

library(survival)
library(survRM2)

library(foreach)
library(doParallel)

# Cluster object
cl = parallel::makeCluster(parallel::detectCores() - 2)
# Register parallel backend
registerDoParallel(cl)

# Set background to be white for all ggplots
theme_set(theme_classic())

Data processing

We read in the data and extract the outcome and covariate variables of interest. Duplicate observations and observations where the treatment arm or variables needed to calculate the KFRE are missing are dropped.

# Read in datasets
# Raw data
raw_data = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/Jane-Vivek/ACCORD HTE Analysis/Data/New data 26SEP2022/accordall.csv", header = TRUE)

# Additional outcomes
other_df = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/Jane-Vivek/ACCORD HTE Analysis/Data/ACCORD_GLY_NEWOUTCOMES-2.csv", header = TRUE)
hypo_df = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/3-Data_Sets-Analysis/3a-Analysis_Data_Sets/csv/hypoglycemiatime1st.csv")
cvd_df = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/3-Data_Sets-Analysis/3a-Analysis_Data_Sets/csv/cvdoutcomes.csv")

# Concomitant medications
meds_df = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/3-Data_Sets-Analysis/3a-Analysis_Data_Sets/csv/concomitantmeds.csv")
# Classification groups for medication
meds_groups_df = read.csv("Copy of Meds_ACCORD.csv")[,1:2]
names(meds_groups_df) = c("med", "group")


# Subset for death outcomes
other_df = other_df %>%
  select(maskid = id, ALLDEATH, alldeath_fu, CVDEATH, cvdeath_fu)
# Subset for 1st assisted hypoglycemic event
hypo_df = hypo_df %>% 
  mutate(hypoglycemia = 1 - censor_any, 
         hypoglycemia_fu = 365 * fuyrs_any) %>% 
  select(maskid = MaskID, hypoglycemia, hypoglycemia_fu)
# Subset for CVD primary outcome
cvd_df = cvd_df %>% 
  mutate(cvd_primary = 1 - censor_po, 
         cvd_primary_fu = 365 * fuyrs_po) %>% 
  select(maskid = MaskID, cvd_primary, cvd_primary_fu)

# Use friendlier variable names
meds_groups_df = meds_groups_df %>% 
  filter(group != "") %>% 
  mutate(group = recode(group, "anti-htn" = "anti_htn", 
                        "chol lowering" = "chol_lowering", 
                        "oral DM" = "oral_DM"))
# Convert data frame of classification groups to list
con_meds = sort(unique(meds_groups_df$group))
meds_groups_list = sapply(con_meds, function(x) {
  meds_groups_df$med[meds_groups_df$group == x]
})
# Only consider baseline meds
meds_df = meds_df %>% filter(Visit == "BLR")
# Create a new variable for each con med classification group
for (med in con_meds) {
  meds_df[[med]] = 
    rowSums(meds_df[,meds_groups_list[[med]]], na.rm = TRUE) > 0
}


# Merge in death outcomes
raw_data = merge(raw_data, other_df, by = "maskid", all.x = TRUE)
# Merge in 1st assisted hypoglycemic event
raw_data = merge(raw_data, hypo_df, by = "maskid", all.x = TRUE)
# Merge in CVD outcomes
raw_data = merge(raw_data, cvd_df, by = "maskid", all.x = TRUE)
# Merge in con meds
raw_data = merge(raw_data %>% select(!diuretic), 
                 meds_df %>% select(maskid = MaskID, all_of(con_meds)), 
                 by = "maskid", all.x = TRUE)

# Create composite kidney outcome
raw_data$neph_composite = pmax(raw_data$Neph2, raw_data$Neph3)
raw_data$neph_composite_fu = pmin(raw_data$Neph2Days, raw_data$Neph3Days)

# Quantitative variables to potentially include in the models 
subset_quant = c("female", "baseline_age", "raceth", 
                 "sbp", "dbp", "pulsepres",
                 "chol", "trig", "vldl", "ldl", "hdl",
                 "alt", "cpk",
                 "fpg", "hba1c", 
                 "ualb", "ucreat", "uacr", 
                 "ckd2021GFR", "screat", 
                 "bmi", 
                 "smokelif", 
                 con_meds)

# List of outcome variables 
outcome_list = list(
  neph_composite = c(status = "neph_composite", time = "neph_composite_fu"), 
  Neph2 = c(status = "Neph2", time = "Neph2Days"), 
  Neph3 = c(status = "Neph3", time = "Neph3Days"), 
  ALLDEATH = c(status = "ALLDEATH", time = "alldeath_fu"), 
  CVDEATH = c(status = "CVDEATH", time = "cvdeath_fu"), 
  hypoglycemia = c(status = "hypoglycemia", time = "hypoglycemia_fu"), 
  cvd_primary = c(status = "cvd_primary", time = "cvd_primary_fu"))
# All outcomes
outcomes = unlist(outcome_list, use.names = FALSE)

# Race should be a factor
raw_data$raceth = as.factor(raw_data$raceth)

# Make relevant exclusions and define new dataset 
df = raw_data %>% 
  select(all_of(c("maskid", "glyarm", "kfrs", outcomes, subset_quant))) %>%
  drop_na(all_of(c("glyarm", "kfrs"))) %>% # Drop those with missing values
  unique() # Drop duplicates

# Include race contrasts
df = data.frame(df, 
                raceth0 = ifelse(df$raceth == 0, 1, 0), 
                model.matrix(~raceth, df)[,-1])
# Read in the longitudinal labs data and calculate eGFR. 
dir = "/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/3-Data_Sets-Analysis/3a-Analysis_Data_Sets/csv/"
activitystatus = read.csv(paste0(dir, "activitystatus.csv"))
otherlabs = read.csv(paste0(dir, "otherlabs.csv"))
accord_key = read.csv(paste0(dir, "accord_key.csv"))

# Calculate eGFR
calc_eGFR = function(AGE, MALE, SCR) {
  if (length(AGE) !=  length(MALE) || 
      length(AGE) !=  length(SCR) || 
      length(MALE) !=  length(SCR)) {
    warning("Input arguments are not the same length.")
  }
  
  kau = 0.7*(1 - MALE) + 0.9*MALE
  alpha = -0.241*(1 - MALE) - 0.302*MALE
  beta = 1.012*(1 - MALE) + 1*MALE
  
  return(142 * pmin(SCR/kau, 1)^alpha * pmax(SCR/kau, 1)^(-1.2) * 
           0.9938^AGE * beta)
}

# Combine otherlabs (contains serum creatinine) with activitystatus (contains 
# FU time) and accord_key (contains age, sex, and arm)
# Calculate eGFR and 6-month spline
# Only include observations with 7 years and for IDs included in the main df
egfr_df = right_join(activitystatus %>% select(MaskID, Visit, days_from_baseline), 
                     otherlabs %>% select(MaskID, Visit, screat, gfr), 
                     by = c("MaskID", "Visit")) %>% 
  left_join(accord_key %>% 
              select(MaskID, female, baseline_age, raceclass, arm, treatment), 
            by = "MaskID") %>% 
  mutate(treat = ifelse(arm %in% c(3, 4, 7, 8), 1, 0), 
         baseline_age = baseline_age + days_from_baseline / 365) %>%
  mutate(egfr = calc_eGFR(baseline_age, female == 0, screat), 
         spline = ifelse(days_from_baseline >= 180, 
                         days_from_baseline - 180, 0)) %>% 
  filter(days_from_baseline <= 365*7, 
         MaskID %in% df$maskid) %>% 
  drop_na()

# Create data frame for time to reduction of eGFR by 40% of baseline
egfr40_df = left_join(
  egfr_df %>% 
    group_by(MaskID) %>% 
    arrange(desc(days_from_baseline)) %>% 
    slice(1) %>% 
    ungroup() %>% 
    select(MaskID, days_from_baseline), 
  egfr_df %>% 
    group_by(MaskID) %>% 
    arrange(days_from_baseline) %>% 
    mutate(egfr40_cumsum = cumsum(egfr <= egfr[days_from_baseline == 0] * 0.6), 
           egfr40_consec_sum = rollsum(egfr <= egfr[days_from_baseline == 0] * 0.6, 
                                       k = 2, na.pad = TRUE,  align = "right")) %>% 
    summarize(
      # 40% eGFR reduction for 2 consecutive measures
      egfr40_2consec = as.numeric(any(egfr40_consec_sum == 2, na.rm = TRUE)),
      egfr40_2consec_fu = first(days_from_baseline[egfr40_consec_sum == 2], 
                                na_rm = TRUE)) %>% 
    ungroup(), 
  by = "MaskID") %>% 
  rename(maskid = MaskID) %>% 
  mutate(egfr40_2consec = replace_na(egfr40_2consec, 0), 
         egfr40_2consec_fu = coalesce(egfr40_2consec_fu, days_from_baseline))

# Merge in event status and time to reduction of eGFR by 40% of baseline
df = left_join(df, egfr40_df, by = "maskid") 

# Create composite out for first occurrence of 40% eGFR reduction or Neph3
# 40% eGFR reduction for 2 consecutive measures
df$egfr40_2consec = pmax(df$egfr40_2consec, df$Neph3)
df$egfr40_2consec_fu = pmin(df$egfr40_2consec_fu, df$Neph3Days)
  
# List of outcome variables 
outcome_list = c(
  outcome_list, 
  list(egfr40_2consec = c(status = "egfr40_2consec", time = "egfr40_2consec_fu")))
# All outcomes
outcomes = unlist(outcome_list, use.names = FALSE)

We define a data frame that standardizes all continuous variables to have mean 0 and standard deviation 1.

# All covariates considered + race split into contrasts
expanded_covariates = c(setdiff(subset_quant, "raceth"), 
                        paste0("raceth", 0:3))
# Binary variables
binary_vars = apply(df[,expanded_covariates], 2, function(x){
  length(unique(x)) == 2 && all(sort(unique(x)) == c(0, 1))
})

# Data where continuous covariates are standardized
# Binary variables are left alone
std_df = sapply(1:length(expanded_covariates), function(i) {
  if (binary_vars[i] == FALSE) {
    return(scale(df[,expanded_covariates][,i]))
  } else {
    return(df[,expanded_covariates][,i])
  }
})
colnames(std_df) = expanded_covariates

Kidney failure risk equation (KFRE)

We calculate the 5-year predicted risk using the KFRE, define KFRE quartiles, and plot the distributions within treatment arms.

# Calculate 5-year KFRE
df$kfre5 = with(df, {
  1 - 0.8996^exp(-0.2201 * (baseline_age/10 - 7.036) + 
                   0.2467 * ((female==0) - 0.5642) - 
                   0.5567 * (ckd2021GFR/5 - 7.222) + 
                   0.4510 * ( log(uacr) - 5.137))
})

# kfre quartile thresholds
kfre_quart_thresh = quantile(df$kfre5, seq(0.25, 0.75, by = 0.25))
# Create groups based on thresholds
df$kfre_quarts = as.numeric(cut(df$kfre5, 
                                breaks = c(-Inf, as.numeric(kfre_quart_thresh), Inf),
                                labels = c(1:(length(kfre_quart_thresh)+1))))

# Histograms of KFRE for treatment and control
df %>% 
  ggplot(aes(x = kfre5, fill = as.factor(glyarm), color = as.factor(glyarm))) + 
  geom_histogram(alpha = 0.3, position="identity") + 
  scale_x_log10(breaks = 10^(-6:0), 
                labels = function(x) paste0(formatC(x*100, format = "fg"), "%")) + 
  geom_vline(xintercept = kfre_quart_thresh, color = "gray60") + 
  geom_text(data = data.frame(kfre_quart_thresh = kfre_quart_thresh), 
             aes(x = kfre_quart_thresh, y = Inf, 
                 label = sprintf("%1.3f%%", 100*kfre_quart_thresh)), 
             inherit.aes = FALSE, color = "grey60", 
            hjust = 1.1, vjust = 1.4, size = 3, angle = 90) + 
  xlab("5-year predicted risk by KFRE (log)") + ylab("") + 
  scale_fill_discrete(name = "", labels = c("0" = "Control", "1" = "Treatment")) + 
  scale_color_discrete(name = "", labels = c("0" = "Control", "1" = "Treatment"))

ggsave("figs_and_tabs/kfre_hist.png", width = 6, height = 4)

# General summary function that returns the mean, standard deviation, median 
# and IQR of x
my_summary = function(x, na.rm = TRUE) {
  c("Mean" = mean(x, na.rm = na.rm), 
    "SD" = sd(x, na.rm = na.rm), 
    "Median" = median(x, na.rm = na.rm), 
    "IQR" = IQR(x, na.rm = na.rm))
}

# Distribution of KFRE for treatment and control
kfre_dist_tab = data.frame(
  do.call(rbind, c(list(my_summary(df$kfre5)), 
                 tapply(df$kfre5, df$glyarm, my_summary))), 
  row.names = c("Overall", "Control", "Treatment")
)
write.csv(kfre_dist_tab %>% round(5), file = "figs_and_tabs/kfre_dist_tab.csv")
kfre_dist_tab %>% round(5)
##              Mean      SD  Median     IQR
## Overall   0.00259 0.01875 0.00014 0.00075
## Control   0.00226 0.01229 0.00013 0.00071
## Treatment 0.00291 0.02346 0.00015 0.00079

These are summary statistics for the scores within each KFRE quartile.

# Summary stats for KFRE within quartile
kfre_quart_tab = data.frame(do.call(rbind, 
        lapply(1:4, function(i) {
          c(Quartile = i, 
            my_summary(df$kfre5[df$kfre_quarts==i]))
        })
))
write.csv(kfre_quart_tab %>% round(5), 
          file = "figs_and_tabs/kfre_quart_tab.csv", row.names = FALSE)
kfre_quart_tab %>% round(5)
##   Quartile    Mean      SD  Median     IQR
## 1        1 0.00002 0.00001 0.00002 0.00001
## 2        2 0.00008 0.00003 0.00007 0.00005
## 3        3 0.00034 0.00017 0.00029 0.00026
## 4        4 0.00991 0.03653 0.00275 0.00575

Number of observations overall and in each arm.

c(table(raw_data$glyarm), Overall = nrow(raw_data))
##       0       1 Overall 
##    5126    5140   10266

Number of observations overall and in each arm, restricting to those with non-missing covariates needed to calculate KFRE.

c(table(df$glyarm), Overall = nrow(df))
##       0       1 Overall 
##    4873    4904    9777

Number of non-missing observations for each outcome.

t(sapply(outcome_list, function(x) {
  my_df = na.omit(df[,c(x["status"], x["time"], "glyarm")])
  c(table(my_df$glyarm), Overall = nrow(my_df))
}))
##                   0    1 Overall
## neph_composite 4237 4219    8456
## Neph2          4237 4219    8456
## Neph3          4865 4895    9760
## ALLDEATH       4873 4904    9777
## CVDEATH        4873 4904    9777
## hypoglycemia   4873 4904    9777
## cvd_primary    4873 4904    9777
## egfr40_2consec 4865 4895    9760

Descriptive statistics

For each quartile and overall, we report the proportion (standard deviation) for all binary variables and the median (IQR) of all continuous variables in the dataset.

# Data dictionary for nice variable names
dat_dict = data.frame(rbind(
  c("female", "Female gender", "(yes/no)"), 
  c("raceth0", "White", "(yes/no)"), 
  c("raceth1", "Black", "(yes/no)"), 
  c("raceth3", "Hispanic", "(yes/no)"), 
  c("raceth2", "Other race/ethnicity", "(yes/no)"), 
  c("baseline_age", "Baseline age", "(years)"), 
  c("bmi", "Body mass index", "(kg/m2)"), 
  c("sbp", "Systolic blood pressure", "(mmHg)"), 
  c("dbp", "Diastolic blood pressure", "(mmHg)"), 
  c("pulsepres", "Pulse pressure", "(mmHg)"), 
  c("chol", "Total cholesterol", "(mg/dL)"), 
  c("trig", "Triglyercides", "(mg/dL)"), 
  c("vldl", "Very low density lipoprotein", "(mg/dL)"), 
  c("ldl", "Low density lipoprotein", "(mg/dL)"), 
  c("hdl", "High density lipoprotein", "(mg/dL)"), 
  c("alt", "Alanine aminotransferase", "(mg/dL)"), 
  c("cpk", "Creatine phosphokinase", "(mg/dL)"), 
  c("fpg", "Fasting plasma glucose", "(mg/dL)"), 
  c("hba1c", "Glycosylated hemoglobin", "(%)"), 
  c("ualb", "Urinary albumin", "(mg/dL)"), 
  c("ucreat", "Urinary creatinine", "(mg/dL)"), 
  c("uacr", "Urinary albumin to creatinine ratio", "(mg/g)"), 
  c("ckd2021GFR", "Estimated glomerular filtration rate", "(mL/min/1.73m2)"), 
  c("screat", "Serum creatinine", "(mg/dL)"),  
  c("anti_htn", "Anti-hypertension medication", "(yes/no)"), 
  c("chol_lowering", "Cholesterol-lowering medication", "(yes/no)"), 
  c("diuretic", "Diuretic medication", "(yes/no)"), 
  c("insulin", "Insulin", "(yes/no)"), 
  c("oral_DM", "Oral diabetes medication", "(yes/no)"), 
  c("neph_composite", "Composite kidney outcome", "(yes/no)"), 
  c("Neph2", "Development of macro-albuminuria", "(yes/no)"), 
  c("Neph3", "Renal failure or ESRD (dialysis) or SCr>3.3", "(yes/no)"), 
  c("ALLDEATH", "All cause death", "(yes/no)"), 
  c("CVDEATH", "Cardiovascular death", "(yes/no)"), 
  c("hypoglycemia", "1st assisted hypoglycemic event", "(yes/no)"), 
  c("egfr40_2consec", "ESRD or sustained (consecutive) eGFR reduction of 40%", "(yes/no)")
))
names(dat_dict) = c("short", "long", "units")

# Center (median or mean) of each variable
center_tab = cbind(
  # Overall
  sapply(c(expanded_covariates, names(outcome_list)), function(var) {
    if (var %in% c("glyarm", names(outcome_list), 
                   names(binary_vars)[binary_vars])) { # Binary variables
      mean(df[,var], na.rm = TRUE)
    } else { # Continuous
      median(df[,var], na.rm = TRUE)
    }
  }), 
  # By quartile
  sapply(1:4, function(i) {
    sapply(c(expanded_covariates, names(outcome_list)), function(var) {
      if (var %in% c("glyarm", names(outcome_list), 
                     names(binary_vars)[binary_vars])) { # Binary variables
        mean(df[df$kfre_quarts==i,var], na.rm = TRUE)
      } else { # Continuous
        median(df[df$kfre_quarts==i,var], na.rm = TRUE)
      }
    })
  })
)

# Spread (IQR or standard deviation) of each variable
spread_tab = cbind(
  # Overall
  sapply(c(expanded_covariates, names(outcome_list)), function(var) {
    if (var %in% c("glyarm", names(outcome_list), 
                   names(binary_vars)[binary_vars])) { # Binary variables
      sd(df[,var], na.rm = TRUE)
    } else { # Continuous
      IQR(df[,var], na.rm = TRUE)
    }
  }), 
  # By quartile
  sapply(1:4, function(i){
    sapply(c(expanded_covariates, names(outcome_list)), function(var) {
      if (var %in% c("glyarm", names(outcome_list), 
                     names(binary_vars)[binary_vars])) { # Binary variables
        sd(df[df$kfre_quarts==i,var], na.rm = TRUE)
      } else { # Continuous
        IQR(df[df$kfre_quarts==i,var], na.rm = TRUE)
      }
    })
  })
)

# Number of non-missing values for each covariate
N_tab = sapply(c(expanded_covariates, names(outcome_list)), function(var) {
  sum(!is.na(df[,var]))
})

# Create Table 1
tab1 = matrix(paste0(round(center_tab, 2), " (", round(spread_tab, 2), ")"), 
              nrow = nrow(center_tab), ncol = 5)
rownames(tab1) = rownames(center_tab)
colnames(tab1) = c("Overall", paste("Quartile", 1:4))
tab1 = cbind(tab1, N = N_tab)

# Use more descriptive variable names
tab1 = merge(dat_dict, tab1, by.x = "short", by.y = 0, sort = FALSE)
rownames(tab1) = paste(tab1$long, tab1$units)
tab1$short = NULL; tab1$long = NULL; tab1$units = NULL
write.csv(tab1, file = "figs_and_tabs/tab1.csv")
tab1
##                                                                      Overall
## Female gender (yes/no)                                           0.38 (0.49)
## White (yes/no)                                                   0.62 (0.48)
## Black (yes/no)                                                   0.19 (0.39)
## Hispanic (yes/no)                                                0.07 (0.26)
## Other race/ethnicity (yes/no)                                    0.11 (0.32)
## Baseline age (years)                                              62.1 (9.5)
## Body mass index (kg/m2)                                         31.83 (7.71)
## Systolic blood pressure (mmHg)                                      136 (22)
## Diastolic blood pressure (mmHg)                                      75 (14)
## Pulse pressure (mmHg)                                                60 (19)
## Total cholesterol (mg/dL)                                           178 (53)
## Triglyercides (mg/dL)                                              156 (122)
## Very low density lipoprotein (mg/dL)                                 31 (25)
## Low density lipoprotein (mg/dL)                                     101 (44)
## High density lipoprotein (mg/dL)                                     40 (14)
## Alanine aminotransferase (mg/dL)                                     24 (14)
## Creatine phosphokinase (mg/dL)                                      106 (92)
## Fasting plasma glucose (mg/dL)                                      168 (65)
## Glycosylated hemoglobin (%)                                        8.1 (1.3)
## Urinary albumin (mg/dL)                                          1.69 (4.48)
## Urinary creatinine (mg/dL)                                      116.9 (77.5)
## Urinary albumin to creatinine ratio (mg/g)                           14 (40)
## Estimated glomerular filtration rate (mL/min/1.73m2)           87.95 (26.47)
## Serum creatinine (mg/dL)                                           0.9 (0.2)
## Anti-hypertension medication (yes/no)                            0.81 (0.39)
## Cholesterol-lowering medication (yes/no)                         0.68 (0.47)
## Diuretic medication (yes/no)                                     0.37 (0.48)
## Insulin (yes/no)                                                 0.35 (0.48)
## Oral diabetes medication (yes/no)                                0.83 (0.37)
## Composite kidney outcome (yes/no)                                 0.1 (0.29)
## Development of macro-albuminuria (yes/no)                        0.07 (0.26)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no)             0.03 (0.17)
## All cause death (yes/no)                                         0.07 (0.26)
## Cardiovascular death (yes/no)                                    0.03 (0.18)
## 1st assisted hypoglycemic event (yes/no)                         0.12 (0.33)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no)   0.11 (0.31)
##                                                                   Quartile 1
## Female gender (yes/no)                                            0.48 (0.5)
## White (yes/no)                                                   0.66 (0.47)
## Black (yes/no)                                                   0.11 (0.31)
## Hispanic (yes/no)                                                0.09 (0.28)
## Other race/ethnicity (yes/no)                                    0.14 (0.35)
## Baseline age (years)                                                58.6 (6)
## Body mass index (kg/m2)                                         32.27 (8.13)
## Systolic blood pressure (mmHg)                                      133 (20)
## Diastolic blood pressure (mmHg)                                      76 (13)
## Pulse pressure (mmHg)                                                56 (17)
## Total cholesterol (mg/dL)                                           182 (54)
## Triglyercides (mg/dL)                                              161 (131)
## Very low density lipoprotein (mg/dL)                                 32 (26)
## Low density lipoprotein (mg/dL)                                     102 (45)
## High density lipoprotein (mg/dL)                                     40 (13)
## Alanine aminotransferase (mg/dL)                                     25 (14)
## Creatine phosphokinase (mg/dL)                                       91 (74)
## Fasting plasma glucose (mg/dL)                                      173 (67)
## Glycosylated hemoglobin (%)                                        8.1 (1.3)
## Urinary albumin (mg/dL)                                           0.9 (1.09)
## Urinary creatinine (mg/dL)                                      112.7 (73.4)
## Urinary albumin to creatinine ratio (mg/g)                             8 (9)
## Estimated glomerular filtration rate (mL/min/1.73m2)           102.37 (6.11)
## Serum creatinine (mg/dL)                                           0.7 (0.2)
## Anti-hypertension medication (yes/no)                            0.74 (0.44)
## Cholesterol-lowering medication (yes/no)                         0.65 (0.48)
## Diuretic medication (yes/no)                                     0.25 (0.44)
## Insulin (yes/no)                                                 0.28 (0.45)
## Oral diabetes medication (yes/no)                                0.85 (0.36)
## Composite kidney outcome (yes/no)                                0.05 (0.21)
## Development of macro-albuminuria (yes/no)                        0.02 (0.13)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no)             0.03 (0.16)
## All cause death (yes/no)                                         0.04 (0.19)
## Cardiovascular death (yes/no)                                    0.01 (0.12)
## 1st assisted hypoglycemic event (yes/no)                         0.08 (0.27)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no)   0.08 (0.27)
##                                                                  Quartile 2
## Female gender (yes/no)                                           0.3 (0.46)
## White (yes/no)                                                  0.67 (0.47)
## Black (yes/no)                                                  0.14 (0.35)
## Hispanic (yes/no)                                               0.07 (0.26)
## Other race/ethnicity (yes/no)                                   0.12 (0.33)
## Baseline age (years)                                               62 (8.6)
## Body mass index (kg/m2)                                        31.81 (7.47)
## Systolic blood pressure (mmHg)                                     135 (21)
## Diastolic blood pressure (mmHg)                                     75 (14)
## Pulse pressure (mmHg)                                               60 (18)
## Total cholesterol (mg/dL)                                          177 (53)
## Triglyercides (mg/dL)                                           156 (121.5)
## Very low density lipoprotein (mg/dL)                                31 (25)
## Low density lipoprotein (mg/dL)                                     99 (43)
## High density lipoprotein (mg/dL)                                    39 (13)
## Alanine aminotransferase (mg/dL)                                    25 (15)
## Creatine phosphokinase (mg/dL)                                     102 (83)
## Fasting plasma glucose (mg/dL)                                     170 (63)
## Glycosylated hemoglobin (%)                                       8.1 (1.3)
## Urinary albumin (mg/dL)                                         1.94 (4.14)
## Urinary creatinine (mg/dL)                                     116.3 (80.1)
## Urinary albumin to creatinine ratio (mg/g)                          17 (40)
## Estimated glomerular filtration rate (mL/min/1.73m2)           95.91 (9.16)
## Serum creatinine (mg/dL)                                          0.8 (0.2)
## Anti-hypertension medication (yes/no)                             0.8 (0.4)
## Cholesterol-lowering medication (yes/no)                        0.68 (0.47)
## Diuretic medication (yes/no)                                    0.32 (0.46)
## Insulin (yes/no)                                                 0.3 (0.46)
## Oral diabetes medication (yes/no)                               0.85 (0.35)
## Composite kidney outcome (yes/no)                               0.09 (0.29)
## Development of macro-albuminuria (yes/no)                       0.07 (0.25)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no)            0.02 (0.15)
## All cause death (yes/no)                                        0.06 (0.24)
## Cardiovascular death (yes/no)                                   0.02 (0.15)
## 1st assisted hypoglycemic event (yes/no)                        0.09 (0.29)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no)  0.11 (0.31)
##                                                                   Quartile 3
## Female gender (yes/no)                                           0.34 (0.47)
## White (yes/no)                                                   0.61 (0.49)
## Black (yes/no)                                                   0.23 (0.42)
## Hispanic (yes/no)                                                0.06 (0.25)
## Other race/ethnicity (yes/no)                                      0.1 (0.3)
## Baseline age (years)                                              63.6 (9.4)
## Body mass index (kg/m2)                                          31.7 (7.64)
## Systolic blood pressure (mmHg)                                      136 (23)
## Diastolic blood pressure (mmHg)                                      75 (15)
## Pulse pressure (mmHg)                                                61 (18)
## Total cholesterol (mg/dL)                                           176 (51)
## Triglyercides (mg/dL)                                              148 (115)
## Very low density lipoprotein (mg/dL)                                 30 (23)
## Low density lipoprotein (mg/dL)                                     100 (42)
## High density lipoprotein (mg/dL)                                     40 (13)
## Alanine aminotransferase (mg/dL)                                  24 (13.12)
## Creatine phosphokinase (mg/dL)                                     113 (102)
## Fasting plasma glucose (mg/dL)                                   167 (63.25)
## Glycosylated hemoglobin (%)                                        8.1 (1.2)
## Urinary albumin (mg/dL)                                          1.85 (5.54)
## Urinary creatinine (mg/dL)                                     121.3 (78.55)
## Urinary albumin to creatinine ratio (mg/g)                         15 (44.5)
## Estimated glomerular filtration rate (mL/min/1.73m2)            81.83 (9.54)
## Serum creatinine (mg/dL)                                             1 (0.2)
## Anti-hypertension medication (yes/no)                            0.83 (0.37)
## Cholesterol-lowering medication (yes/no)                          0.7 (0.46)
## Diuretic medication (yes/no)                                     0.39 (0.49)
## Insulin (yes/no)                                                 0.37 (0.48)
## Oral diabetes medication (yes/no)                                0.84 (0.37)
## Composite kidney outcome (yes/no)                                0.09 (0.29)
## Development of macro-albuminuria (yes/no)                        0.07 (0.26)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no)             0.02 (0.15)
## All cause death (yes/no)                                         0.07 (0.25)
## Cardiovascular death (yes/no)                                    0.03 (0.18)
## 1st assisted hypoglycemic event (yes/no)                         0.13 (0.34)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no)   0.11 (0.32)
##                                                                    Quartile 4
## Female gender (yes/no)                                            0.39 (0.49)
## White (yes/no)                                                     0.57 (0.5)
## Black (yes/no)                                                    0.29 (0.45)
## Hispanic (yes/no)                                                 0.06 (0.24)
## Other race/ethnicity (yes/no)                                     0.09 (0.28)
## Baseline age (years)                                             65.15 (10.3)
## Body mass index (kg/m2)                                           31.54 (7.6)
## Systolic blood pressure (mmHg)                                       138 (23)
## Diastolic blood pressure (mmHg)                                       73 (16)
## Pulse pressure (mmHg)                                                 65 (21)
## Total cholesterol (mg/dL)                                            178 (52)
## Triglyercides (mg/dL)                                               158 (121)
## Very low density lipoprotein (mg/dL)                                  32 (24)
## Low density lipoprotein (mg/dL)                                    101 (44.5)
## High density lipoprotein (mg/dL)                                      40 (15)
## Alanine aminotransferase (mg/dL)                                      22 (12)
## Creatine phosphokinase (mg/dL)                                      120 (113)
## Fasting plasma glucose (mg/dL)                                       163 (67)
## Glycosylated hemoglobin (%)                                         8.2 (1.3)
## Urinary albumin (mg/dL)                                          3.99 (15.82)
## Urinary creatinine (mg/dL)                                     116.75 (76.23)
## Urinary albumin to creatinine ratio (mg/g)                           35 (153)
## Estimated glomerular filtration rate (mL/min/1.73m2)             63.3 (13.24)
## Serum creatinine (mg/dL)                                            1.2 (0.3)
## Anti-hypertension medication (yes/no)                             0.87 (0.33)
## Cholesterol-lowering medication (yes/no)                          0.71 (0.46)
## Diuretic medication (yes/no)                                       0.51 (0.5)
## Insulin (yes/no)                                                   0.45 (0.5)
## Oral diabetes medication (yes/no)                                 0.78 (0.41)
## Composite kidney outcome (yes/no)                                 0.17 (0.37)
## Development of macro-albuminuria (yes/no)                         0.15 (0.36)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no)               0.04 (0.2)
## All cause death (yes/no)                                          0.12 (0.32)
## Cardiovascular death (yes/no)                                     0.06 (0.24)
## 1st assisted hypoglycemic event (yes/no)                          0.18 (0.39)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no)    0.12 (0.33)
##                                                                   N
## Female gender (yes/no)                                         9777
## White (yes/no)                                                 9777
## Black (yes/no)                                                 9777
## Hispanic (yes/no)                                              9777
## Other race/ethnicity (yes/no)                                  9777
## Baseline age (years)                                           9777
## Body mass index (kg/m2)                                        9771
## Systolic blood pressure (mmHg)                                 9733
## Diastolic blood pressure (mmHg)                                9733
## Pulse pressure (mmHg)                                          9733
## Total cholesterol (mg/dL)                                      9775
## Triglyercides (mg/dL)                                          9775
## Very low density lipoprotein (mg/dL)                           9775
## Low density lipoprotein (mg/dL)                                9775
## High density lipoprotein (mg/dL)                               9776
## Alanine aminotransferase (mg/dL)                               9777
## Creatine phosphokinase (mg/dL)                                 9776
## Fasting plasma glucose (mg/dL)                                 9777
## Glycosylated hemoglobin (%)                                    9763
## Urinary albumin (mg/dL)                                        9777
## Urinary creatinine (mg/dL)                                     9777
## Urinary albumin to creatinine ratio (mg/g)                     9777
## Estimated glomerular filtration rate (mL/min/1.73m2)           9777
## Serum creatinine (mg/dL)                                       9777
## Anti-hypertension medication (yes/no)                          9777
## Cholesterol-lowering medication (yes/no)                       9777
## Diuretic medication (yes/no)                                   9777
## Insulin (yes/no)                                               9777
## Oral diabetes medication (yes/no)                              9777
## Composite kidney outcome (yes/no)                              8456
## Development of macro-albuminuria (yes/no)                      8456
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no)           9777
## All cause death (yes/no)                                       9777
## Cardiovascular death (yes/no)                                  9777
## 1st assisted hypoglycemic event (yes/no)                       9777
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no) 9777

Balance between treatment and control arm within quartile

To check the balance between the treatment and control groups, we take the difference in standardized means for each covariate, within each treatment arm and quartile.

# Mean of each covariate across individuals within each quartile and treatment 
# group, for kfirs
# Standardized means were used for continuous covariates
std_means_balance_kfre = lapply(expanded_covariates, function(var) {
  tapply(std_df[,var], list(df$kfre_quarts, df$glyarm), 
         mean, na.rm = TRUE)
})
names(std_means_balance_kfre) = expanded_covariates

# Split into treatment and control
# Treatment
std_means_trt_kfre = sapply(expanded_covariates, function(var){
  std_means_balance_kfre[[var]][,"1"]
}) %>% t()

# Control
std_means_control_kfre = sapply(expanded_covariates, function(var){
  std_means_balance_kfre[[var]][,"0"]
}) %>% t()

# Difference between treatment and control
std_means_diff_kfre = 
  merge(dat_dict, 
        sapply(expanded_covariates, function(var){
          std_means_balance_kfre[[var]][,"1"] - 
            std_means_balance_kfre[[var]][,"0"]
          }) %>% t(), 
        by.x = "short", by.y = 0, sort = FALSE)
rownames(std_means_diff_kfre) = std_means_diff_kfre$long
std_means_diff_kfre = std_means_diff_kfre[,-c(1:3)]

We plot the mean difference between treatment and control for each variable, within each quartile.

data.frame(est = gather(data.frame(std_means_diff_kfre))$value, 
           Quartile = rep(as.factor(1:4), each = nrow(std_means_diff_kfre))) %>% 
  mutate(y = rep(1:nrow(std_means_diff_kfre), 4) + 
           rep(seq(-0.2, 0.2, length = 4), each = nrow(std_means_diff_kfre))) %>% 
  ggplot(aes(x = as.numeric(est), y = y, color = Quartile)) +
  geom_point() + 
  scale_color_manual(values = c("navy", "dodgerblue", "mediumpurple","firebrick")) + 
  scale_y_continuous(breaks = 1:nrow(std_means_diff_kfre), 
                     labels = rownames(std_means_diff_kfre), 
                     expand = c(0.02, 0.02)) + 
  xlab("Difference in standardized mean between treatment and control") + ylab("") + 
  ggtitle("Imbalance by KFRE quartile")

ggsave("figs_and_tabs/kfre_imbalance.png", width = 6, height = 5)

write.csv(std_means_diff_kfre %>% round(3), 
          file = "figs_and_tabs/std_means_diff_kfre.csv")
std_means_diff_kfre %>% round(3)
##                                           1      2      3      4
## Female gender                        -0.019  0.006  0.027  0.010
## White                                -0.010  0.010  0.011 -0.001
## Black                                 0.016  0.002 -0.006  0.004
## Hispanic                              0.009 -0.013 -0.013 -0.002
## Other race/ethnicity                 -0.016  0.002  0.008 -0.001
## Baseline age                         -0.022  0.036 -0.059  0.005
## Body mass index                       0.003 -0.057  0.021 -0.006
## Systolic blood pressure              -0.049  0.034 -0.012 -0.028
## Diastolic blood pressure             -0.025 -0.035  0.005 -0.023
## Pulse pressure                       -0.040  0.066 -0.017 -0.016
## Total cholesterol                     0.066 -0.041  0.005 -0.004
## Triglyercides                         0.021 -0.011  0.051  0.005
## Very low density lipoprotein          0.030 -0.009  0.050 -0.011
## Low density lipoprotein               0.064 -0.044 -0.020 -0.001
## High density lipoprotein             -0.014  0.000 -0.027  0.011
## Alanine aminotransferase              0.021 -0.036  0.009 -0.012
## Creatine phosphokinase                0.010 -0.032 -0.081  0.048
## Fasting plasma glucose                0.083  0.004 -0.064 -0.086
## Glycosylated hemoglobin               0.089 -0.024 -0.076 -0.056
## Urinary albumin                      -0.001 -0.003 -0.022 -0.021
## Urinary creatinine                   -0.064 -0.013 -0.014  0.098
## Urinary albumin to creatinine ratio   0.005 -0.005 -0.006  0.037
## Estimated glomerular filtration rate  0.007  0.003  0.015 -0.026
## Serum creatinine                      0.008 -0.015 -0.025  0.039
## Anti-hypertension medication         -0.031 -0.018  0.006 -0.014
## Cholesterol-lowering medication      -0.027  0.002  0.015 -0.006
## Diuretic medication                  -0.010 -0.006 -0.028  0.021
## Insulin                              -0.017 -0.006 -0.029 -0.021
## Oral diabetes medication             -0.009 -0.009  0.019  0.010

No covariates were imbalanced (defined as the mean difference between treatment and control being greater than 0.1 for at least one quartile), so we will not run an adjusted analysis.

# Which variables are imbalanced by more than 0.05 standard units?
is_imbalanced = apply(std_means_diff_kfre, 1, function(x){
  any(abs(x) > 0.1)
})

imbalanced_vars = dat_dict$short[which(dat_dict$long %in%
                                         names(is_imbalanced)[is_imbalanced])]
imbalanced_vars
## character(0)

Treatment effects

We estimate the 7-year RMST differences for each subgroup defined by the KFRE quartiles.

# Horizon
horizon = 365*7

# RMST difference overall
overall_rmst_fits = lapply(outcome_list, function(x){
  df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
      if (x["status"] == "egfr40_2consec") {
        horizon = 2525 
      } else {
        horizon = 365*7
      }
  rmst2(df_sub[,x["time"]], 
        df_sub[,x["status"]], 
        df_sub[,"glyarm"], 
        tau = horizon)
})
names(overall_rmst_fits) = names(outcome_list)
overall_rmst = sapply(names(outcome_list), function(outcome){
  overall_rmst_fits[[outcome]]$unadjusted.result[1,1]
})

# RMST difference by quartile
kfre_quart_rmst_fits = lapply(outcome_list, function(x){
  df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
  lapply(1:4, function(i){
    rmst2(df_sub[df_sub$kfre_quarts==i,x["time"]], 
          df_sub[df_sub$kfre_quarts==i,x["status"]], 
          df_sub[df_sub$kfre_quarts==i,"glyarm"], 
          tau = horizon)
  })
})
names(kfre_quart_rmst_fits) = names(outcome_list)
kfre_quart_rmst = sapply(names(outcome_list), function(outcome){
  sapply(1:4, function(i){
    kfre_quart_rmst_fits[[outcome]][[i]]$unadjusted.result[1,1]
  })
})

We use 1000 bootstrap samples to obtain confidence intervals for the RMST differences. Tables of the RMST differences and normalized RMST differences (where the overall RMST difference is subtracted from each quartile’s estimate) with 95% bootstrap CIs are shown below.

# Number of bootstraps/permutations
B = 1000
boot_kfre_rmst = foreach::foreach(
  b = 1:B, 
  .packages = c("survRM2")
  ) %dopar% {
    # Set seed
    set.seed(b)
    # Create bootstrap sample
    boot_idx = sample(nrow(df), replace= TRUE)
    df_boot = df[boot_idx,]
    
    # Calculate overall RMST difference
    overall_rmst = sapply(outcome_list, function(x){
      df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
      rmst2(df_sub[,x["time"]], 
            df_sub[,x["status"]], 
            df_sub[,"glyarm"], 
            tau = horizon)$unadjusted.result[1,1]
      })
    
    # Calculate RMST difference for each quartile
    kfre_quart_rmst = sapply(outcome_list, function(x){
      df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
      if (x["status"] == "egfr40_2consec") {
        horizon = 2525 
      } else {
        horizon = 365*7
      }
      sapply(1:4, function(i){
        rmst2(df_sub[df_sub$kfre_quarts==i,x["time"]], 
              df_sub[df_sub$kfre_quarts==i,x["status"]], 
              df_sub[df_sub$kfre_quarts==i,"glyarm"], 
              tau = horizon)$unadjusted.result[1,1]
      })
    })
    return(list(overall_rmst = overall_rmst, 
                kfre_quart_rmst = kfre_quart_rmst))
  }

save(boot_kfre_rmst, file = "boot_kfre_rmst.rData")
# Overall
kfre_overall_rmst_boot_ci = apply(
  sapply(boot_kfre_rmst, function(x){
    x$overall_rmst
  }), 1, quantile, c(0.025, 0.975))

# Not normalized
kfre_quart_rmst_boot_ci = lapply(1:4, function(i){
  dat = sapply(boot_kfre_rmst, function(x){
    x$kfre_quart_rmst[i,]
  })
  apply(dat, 1, quantile, c(0.025, 0.975))
})

# Make tables
kfre_quart_rmst_boot_tabs = lapply(names(outcome_list), function(outcome){
  out = rbind(
    data.frame(Quartile = "Overall", 
               "Est" = overall_rmst[outcome], 
               t(kfre_overall_rmst_boot_ci[,outcome])), 
    data.frame(Quartile = 1:4, 
               "Est" = kfre_quart_rmst[,outcome], 
               t(sapply(kfre_quart_rmst_boot_ci, function(x){x[,outcome]}))))
  rownames(out) = NULL
  return(out)
})
names(kfre_quart_rmst_boot_tabs) = names(outcome_list)

# Format table for printing
rmst_tab = data.frame(
  Outcome = rep(names(kfre_quart_rmst_boot_tabs), 
                each = nrow(kfre_quart_rmst_boot_tabs[[1]])), 
  do.call(rbind, kfre_quart_rmst_boot_tabs), 
  row.names = NULL
)
rmst_tab = merge(dat_dict[,c("short", "long")], rmst_tab, 
                 by.x = "short", by.y = "Outcome", all.y = TRUE, 
                 sort = FALSE)[,-1]
names(rmst_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
rmst_tab[,c("Est", "2.5%", "97.5%")] = 
  round(rmst_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(rmst_tab, 
          file = "figs_and_tabs/rmst_tab.csv", row.names = FALSE)
rmst_tab
##                                                  Outcome Quartile     Est
## 1                               Composite kidney outcome        1   10.14
## 2                               Composite kidney outcome        2   40.07
## 3                               Composite kidney outcome        3   50.08
## 4                               Composite kidney outcome  Overall   48.40
## 5                               Composite kidney outcome        4  114.77
## 6                       Development of macro-albuminuria  Overall   42.49
## 7                       Development of macro-albuminuria        1   10.99
## 8                       Development of macro-albuminuria        2   31.32
## 9                       Development of macro-albuminuria        4  116.49
## 10                      Development of macro-albuminuria        3   33.46
## 11           Renal failure or ESRD (dialysis) or SCr>3.3  Overall    5.43
## 12           Renal failure or ESRD (dialysis) or SCr>3.3        1   -4.10
## 13           Renal failure or ESRD (dialysis) or SCr>3.3        2    7.53
## 14           Renal failure or ESRD (dialysis) or SCr>3.3        4    0.99
## 15           Renal failure or ESRD (dialysis) or SCr>3.3        3   18.22
## 16                                       All cause death  Overall  -23.64
## 17                                       All cause death        1  -21.03
## 18                                       All cause death        2  -19.53
## 19                                       All cause death        3   11.06
## 20                                       All cause death        4  -56.70
## 21                                  Cardiovascular death  Overall  -15.43
## 22                                  Cardiovascular death        2  -16.46
## 23                                  Cardiovascular death        3   11.88
## 24                                  Cardiovascular death        4  -35.01
## 25                                  Cardiovascular death        1  -17.82
## 26                       1st assisted hypoglycemic event  Overall -232.97
## 27                       1st assisted hypoglycemic event        2 -216.06
## 28                       1st assisted hypoglycemic event        3 -247.08
## 29                       1st assisted hypoglycemic event        4 -288.38
## 30                       1st assisted hypoglycemic event        1 -171.59
## 31 ESRD or sustained (consecutive) eGFR reduction of 40%  Overall   -1.86
## 32 ESRD or sustained (consecutive) eGFR reduction of 40%        1   19.97
## 33 ESRD or sustained (consecutive) eGFR reduction of 40%        2   -4.53
## 34 ESRD or sustained (consecutive) eGFR reduction of 40%        3   12.84
## 35 ESRD or sustained (consecutive) eGFR reduction of 40%        4  -29.32
## 36                                                  <NA>  Overall   13.46
## 37                                                  <NA>        1   -9.14
## 38                                                  <NA>        2    6.54
## 39                                                  <NA>        3   41.38
## 40                                                  <NA>        4   23.27
##       2.5%   97.5%
## 1   -23.90   43.13
## 2    -4.65   84.70
## 3     7.95   92.69
## 4    25.34   69.60
## 5    58.14  176.41
## 6    22.52   61.31
## 7    -8.12   30.61
## 8    -4.86   68.63
## 9    62.97  173.51
## 10   -6.22   71.98
## 11   -7.02   16.95
## 12  -30.36   19.86
## 13  -14.43   31.42
## 14  -27.60   27.66
## 15   -3.14   39.30
## 16  -42.20   -6.64
## 17  -47.70    6.66
## 18  -53.51   12.56
## 19  -20.23   42.14
## 20 -100.19  -17.50
## 21  -28.26   -3.57
## 22  -40.91    4.34
## 23  -12.68   39.15
## 24  -65.91   -5.41
## 25  -37.08    0.40
## 26 -259.47 -208.03
## 27 -263.96 -168.89
## 28 -302.38 -192.26
## 29 -346.09 -227.82
## 30 -217.89 -127.53
## 31  -25.85   19.87
## 32  -25.65   66.44
## 33  -52.13   41.24
## 34  -30.46   56.44
## 35  -75.96   16.64
## 36  -10.10   34.64
## 37  -42.10   26.45
## 38  -37.73   48.36
## 39    0.72   86.08
## 40  -33.49   77.24
# Normalized (subtract overall RMST difference)
kfre_quart_rmst_norm_boot_ci = lapply(1:4, function(i){
  dat = sapply(boot_kfre_rmst, function(x){
    x$kfre_quart_rmst[i,] - x$overall_rmst
  })
  apply(dat, 1, quantile, c(0.025, 0.975))
})

# Make tables
kfre_quart_rmst_boot_norm_tabs = lapply(names(outcome_list), function(outcome){
  data.frame(Quartile = 1:4, 
               "Est" = kfre_quart_rmst[,outcome] - overall_rmst[outcome], 
               t(sapply(kfre_quart_rmst_norm_boot_ci, function(x){x[,outcome]})))
})
names(kfre_quart_rmst_boot_norm_tabs) = names(outcome_list)

# Format table for printing
rmst_norm_tab = data.frame(
  Outcome = rep(names(kfre_quart_rmst_boot_norm_tabs), 
                each = nrow(kfre_quart_rmst_boot_norm_tabs[[1]])), 
  do.call(rbind, kfre_quart_rmst_boot_norm_tabs), 
  row.names = NULL
)
rmst_norm_tab = merge(dat_dict[,c("short", "long")], rmst_norm_tab, 
                      by.x = "short", by.y = "Outcome", all.y = TRUE, 
                      sort = FALSE)[,-1]
names(rmst_norm_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
rmst_norm_tab[,c("Est", "2.5%", "97.5%")] = 
  round(rmst_norm_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(rmst_norm_tab, 
          file = "figs_and_tabs/rmst_norm_tab.csv", row.names = FALSE)
rmst_norm_tab
##                                                  Outcome Quartile    Est
## 1                               Composite kidney outcome        1 -38.25
## 2                               Composite kidney outcome        2  -8.33
## 3                               Composite kidney outcome        3   1.68
## 4                               Composite kidney outcome        4  66.37
## 5                       Development of macro-albuminuria        1 -31.50
## 6                       Development of macro-albuminuria        2 -11.17
## 7                       Development of macro-albuminuria        3  -9.02
## 8                       Development of macro-albuminuria        4  74.00
## 9            Renal failure or ESRD (dialysis) or SCr>3.3        1  -9.53
## 10           Renal failure or ESRD (dialysis) or SCr>3.3        2   2.09
## 11           Renal failure or ESRD (dialysis) or SCr>3.3        3  12.79
## 12           Renal failure or ESRD (dialysis) or SCr>3.3        4  -4.45
## 13                                       All cause death        1   2.61
## 14                                       All cause death        2   4.11
## 15                                       All cause death        3  34.70
## 16                                       All cause death        4 -33.06
## 17                                  Cardiovascular death        1  -2.39
## 18                                  Cardiovascular death        2  -1.03
## 19                                  Cardiovascular death        3  27.31
## 20                                  Cardiovascular death        4 -19.58
## 21                       1st assisted hypoglycemic event        1  61.38
## 22                       1st assisted hypoglycemic event        2  16.92
## 23                       1st assisted hypoglycemic event        3 -14.11
## 24                       1st assisted hypoglycemic event        4 -55.41
## 25 ESRD or sustained (consecutive) eGFR reduction of 40%        1  21.82
## 26 ESRD or sustained (consecutive) eGFR reduction of 40%        2  -2.67
## 27 ESRD or sustained (consecutive) eGFR reduction of 40%        3  14.70
## 28 ESRD or sustained (consecutive) eGFR reduction of 40%        4 -27.46
## 29                                                  <NA>        1 -22.61
## 30                                                  <NA>        2  -6.93
## 31                                                  <NA>        3  27.92
## 32                                                  <NA>        4   9.81
##       2.5%  97.5%
## 1   -70.97  -6.26
## 2   -46.88  31.06
## 3   -35.00  36.96
## 4    19.51 118.09
## 5   -55.25  -7.48
## 6   -46.01  21.76
## 7   -42.38  23.75
## 8    28.96 122.59
## 9   -32.56  13.00
## 10  -17.36  24.02
## 11   -7.34  33.15
## 12  -27.63  18.20
## 13  -21.95  30.27
## 14  -27.23  32.29
## 15    7.29  63.81
## 16  -66.43  -0.53
## 17  -20.93  17.13
## 18  -23.66  19.05
## 19    5.59  50.96
## 20  -44.91   4.62
## 21   18.55 103.79
## 22  -28.47  57.71
## 23  -60.68  33.34
## 24 -104.52  -5.58
## 25  -21.43  60.40
## 26  -41.31  38.06
## 27  -24.19  55.58
## 28  -68.53  12.52
## 29  -53.37   9.51
## 30  -44.85  30.90
## 31   -8.48  64.56
## 32  -34.86  54.22

These are plots of the RMST differences by outcome, with a horizontal reference line drawn at the overall/ATE point estimate. Shading denotes the ATE CI.

for (outcome in names(kfre_quart_rmst_boot_tabs)) {
  ate_tab = kfre_quart_rmst_boot_tabs[[outcome]][1,]
  p1 = kfre_quart_rmst_boot_tabs[[outcome]][-1,] %>% 
    mutate(x = 1:4) %>% 
    ggplot(aes(x = x, y = as.numeric(Est))) +
    geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) + 
    geom_hline(yintercept = ate_tab[,"Est"], color = "grey") + 
    annotate("rect", xmin = -Inf, xmax = Inf, 
             ymin = ate_tab[,"X2.5."], ymax = ate_tab[,"X97.5."], 
             alpha = 0.2) + 
    scale_x_continuous(breaks = 1:4, labels = 1:4) + 
    theme(legend.position = 'bottom', legend.box = 'vertical') + 
    xlab("Quartile") + ylab("RMST difference") +
    ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
  print(p1)
  ggsave(paste0("figs_and_tabs/rmst_", outcome, ".png"), 
         width = 6, height = 4)
}

Normalized RMST differences.

for (outcome in names(kfre_quart_rmst_boot_norm_tabs)) {
  p1 = kfre_quart_rmst_boot_norm_tabs[[outcome]] %>% 
    mutate(x = 1:4) %>% 
    ggplot(aes(x = x, y = as.numeric(Est))) +
    geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) + 
    scale_x_continuous(breaks = 1:4, labels = 1:4) + 
    theme(legend.position = 'bottom', legend.box = 'vertical') + 
    xlab("Quartile") + ylab("Normalized RMST difference") + 
    geom_hline(yintercept = 0, color = "grey") + 
    ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
  print(p1)
  ggsave(paste0("figs_and_tabs/rmst_norm_", outcome, ".png"), 
         width = 6, height = 4)
}

Cox PH models and cumulative incidence curves

For each outcome, we fit a Cox PH model to treatment for the entire dataset and each group defined by the KFRE quartiles. Estimates of the treatment hazard ratios are shown below.

# Fit Cox PH model for each outcome to glyarm
coxph_fits = lapply(names(outcome_list), function(outcome){
  # Extract data of interest
  x = outcome_list[[outcome]]
  df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm")])
  
  # Formula for survfit object
  form = formula(paste0("Surv(", x["time"], ", " , 
                        x["status"], ") ~ glyarm"))
  # Fit survfit object
  fit = coxph(form, data = df_sub)
  # Manually insert formula
  fit$call$formula = form
  
  return(fit)
})
names(coxph_fits) = names(outcome_list)

# Fit Cox PH models within quartile
coxph_fits_by_quart = lapply(names(outcome_list), function(outcome){
  # Extract data of interest
  x = outcome_list[[outcome]]
  df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
  
  # Formula for survfit object
  form = formula(paste0("Surv(", x["time"], ", " , 
                        x["status"], ") ~ glyarm"))
  
  all_fits = lapply(1:4, function(i){
    # Subset data to quartile
    df_sub = df_sub[df_sub$kfre_quarts==i, ]
    
    # Fit survfit object
    fit = coxph(form, data = df_sub)
    # Manually insert formula
    fit$call$formula = form
    
    return(fit)
  })
  
  return(all_fits)
})
names(coxph_fits_by_quart) = names(outcome_list)

# Make tables
coxph_coeff_tabs = lapply(names(outcome_list), function(outcome){
  out = rbind(
    data.frame(Quartile = "Overall", 
               t(exp(c(coef(coxph_fits[[outcome]]), confint(coxph_fits[[outcome]]))))
               ), 
    data.frame(Quartile = 1:4, 
               do.call(rbind, 
                       lapply(coxph_fits_by_quart[[outcome]], function(x){
                         exp(c(coef(x), confint(x)))
                         }))
               ))
  rownames(out) = NULL
  colnames(out) = c("Quartile", "Est", "X2.5.", "X97.5.")
  return(out)
})
names(coxph_coeff_tabs) = names(outcome_list)

# Format table for printing
coeff_tab = data.frame(
  Outcome = rep(names(coxph_coeff_tabs), 
                each = nrow(coxph_coeff_tabs[[1]])), 
  do.call(rbind, coxph_coeff_tabs), 
  row.names = NULL
)
coeff_tab = merge(dat_dict[,c("short", "long")], coeff_tab, 
                  by.x = "short", by.y = "Outcome", all.y = TRUE, 
                  sort = FALSE)[,-1]
names(coeff_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
coeff_tab[,c("Est", "2.5%", "97.5%")] = 
  round(coeff_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(coeff_tab, 
          file = "figs_and_tabs/coeff_tab.csv", row.names = FALSE)
coeff_tab
##                                                  Outcome Quartile  Est 2.5%
## 1                               Composite kidney outcome        1 0.93 0.63
## 2                               Composite kidney outcome        2 0.78 0.59
## 3                               Composite kidney outcome        3 0.73 0.54
## 4                               Composite kidney outcome  Overall 0.75 0.65
## 5                               Composite kidney outcome        4 0.66 0.53
## 6                       Development of macro-albuminuria  Overall 0.72 0.61
## 7                       Development of macro-albuminuria        1 0.73 0.39
## 8                       Development of macro-albuminuria        2 0.77 0.56
## 9                       Development of macro-albuminuria        4 0.62 0.49
## 10                      Development of macro-albuminuria        3 0.77 0.56
## 11           Renal failure or ESRD (dialysis) or SCr>3.3  Overall 0.93 0.73
## 12           Renal failure or ESRD (dialysis) or SCr>3.3        1 1.09 0.68
## 13           Renal failure or ESRD (dialysis) or SCr>3.3        2 0.89 0.52
## 14           Renal failure or ESRD (dialysis) or SCr>3.3        4 1.03 0.69
## 15           Renal failure or ESRD (dialysis) or SCr>3.3        3 0.65 0.38
## 16                                       All cause death  Overall 1.20 1.04
## 17                                       All cause death        1 1.37 0.91
## 18                                       All cause death        2 1.15 0.84
## 19                                       All cause death        3 0.95 0.70
## 20                                       All cause death        4 1.30 1.03
## 21                                  Cardiovascular death  Overall 1.30 1.05
## 22                                  Cardiovascular death        2 1.48 0.89
## 23                                  Cardiovascular death        3 0.85 0.55
## 24                                  Cardiovascular death        4 1.40 1.01
## 25                                  Cardiovascular death        1 1.75 0.88
## 26                       1st assisted hypoglycemic event  Overall 2.91 2.56
## 27                       1st assisted hypoglycemic event        2 3.90 2.85
## 28                       1st assisted hypoglycemic event        3 2.99 2.32
## 29                       1st assisted hypoglycemic event        4 2.41 1.97
## 30                       1st assisted hypoglycemic event        1 3.05 2.22
## 31 ESRD or sustained (consecutive) eGFR reduction of 40%  Overall 1.02 0.90
## 32 ESRD or sustained (consecutive) eGFR reduction of 40%        1 0.92 0.70
## 33 ESRD or sustained (consecutive) eGFR reduction of 40%        2 1.03 0.81
## 34 ESRD or sustained (consecutive) eGFR reduction of 40%        3 0.92 0.73
## 35 ESRD or sustained (consecutive) eGFR reduction of 40%        4 1.16 0.93
## 36                                                  <NA>  Overall 0.91 0.81
## 37                                                  <NA>        1 1.08 0.77
## 38                                                  <NA>        2 0.94 0.73
## 39                                                  <NA>        3 0.80 0.62
## 40                                                  <NA>        4 0.89 0.73
##    97.5%
## 1   1.37
## 2   1.04
## 3   0.97
## 4   0.86
## 5   0.83
## 6   0.84
## 7   1.39
## 8   1.06
## 9   0.79
## 10  1.07
## 11  1.18
## 12  1.76
## 13  1.51
## 14  1.53
## 15  1.10
## 16  1.40
## 17  2.07
## 18  1.58
## 19  1.29
## 20  1.64
## 21  1.62
## 22  2.48
## 23  1.31
## 24  1.95
## 25  3.47
## 26  3.30
## 27  5.35
## 28  3.84
## 29  2.94
## 30  4.17
## 31  1.15
## 32  1.22
## 33  1.31
## 34  1.17
## 35  1.46
## 36  1.03
## 37  1.51
## 38  1.23
## 39  1.01
## 40  1.09

We plot cumulative incidence curves and risk tables for each outcome.

# Fit Kaplan-Meier survival curves for each outcome to glyarm
km_fits = lapply(names(outcome_list), function(outcome){
  # Extract data of interest
  x = outcome_list[[outcome]]
  df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
  
  # Formula for survfit object
  form = formula(paste0("Surv(", x["time"], ", " , 
                        x["status"], ") ~ glyarm"))
  # Fit survfit object
  fit = survfit(form, data = df_sub)
  # Manually insert formula
  fit$call$formula = form
  
  return(fit)
})
names(km_fits) = names(outcome_list)

# Helper function for formatting labels of the form est (lo, hi)
est_with_ci_text = function(x, est = "Est", lo = "X2.5.", hi = "X97.5.", 
                            digits = 2) {
  x = round(x[c(est, lo, hi)], 2)
  paste0(x[est], " (", x[lo], ", ", x[hi], ")")
}

invisible(lapply(names(outcome_list), function(outcome){
  # Extract data of interest
  x = outcome_list[[outcome]]
  df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
  
  # Plot cumulative incidence curve with risk table
  if (outcome == "egfr40_2consec") {
    outcome_title = "ESRD or sustained (consecutive) \neGFR reduction of 40%"
  } else {
    outcome_title = dat_dict$long[which(dat_dict$short == outcome)]
  }
  p1 = ggsurvplot(km_fits[[outcome]], data = df_sub, censor = FALSE, 
                  title = outcome_title, 
                  conf.int = TRUE, risk.table = TRUE, fun = "event", 
                  legend.labs = c("Control", "Treatment"))
  # Create annotation for hazard ratio and 7-year RMST
  p1$plot = p1$plot + 
    annotate("label", x = 0, y = Inf, vjust = 1, hjust = 0, 
             size = 4, fontface = 2, 
             label = paste0("HR = ", 
                            est_with_ci_text(coxph_coeff_tabs[[outcome]][1,-1]), 
                            "\n7-year RMST difference\n   = ", 
                            est_with_ci_text(kfre_quart_rmst_boot_tabs[[outcome]][1,-1])))
  
  png(paste0("figs_and_tabs/cum_inc_", outcome, ".png"), 
      res = 100, width = 600, height = 600)
  print(p1)
  dev.off()
  print(p1)
}))

We plot cumulative incidence curves for each outcome, faceted by quartile.

# Fit Kaplan-Meier survival curves within quartile
km_fits_by_quart = lapply(names(outcome_list), function(outcome){
  # Extract data of interest
  x = outcome_list[[outcome]]
  df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
  
  # Formula for survfit object
  form = formula(paste0("Surv(", x["time"], ", " , 
                        x["status"], ") ~ glyarm"))
  
  all_fits = lapply(1:4, function(i){
    # Subset data to quartile
    df_sub = df_sub[df_sub$kfre_quarts==i, ]
    
    # Fit survfit object
    fit = survfit(form, data = df_sub)
    # Manually insert formula
    fit$call$formula = form
    
    return(fit)
  })
  
  return(all_fits)
})
names(km_fits_by_quart) = names(outcome_list)

invisible(lapply(names(outcome_list), function(outcome){
  # Extract data of interest
  x = outcome_list[[outcome]]
  df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
  
  # Formula for survfit object
  form = formula(paste0("Surv(", x["time"], ", " , 
                        x["status"], ") ~ glyarm"))
  
  all_plots = lapply(1:4, function(i){
    # Subset data to quartile
    df_sub = df_sub[df_sub$kfre_quarts==i, ]
    
    # Plot cumulative incidence curve with risk table
    if (outcome == "Neph2") {
      outcome_title = "Development of \nmacro-albuminuria"
    } else if (outcome == "Neph3") {
      outcome_title = "Renal failure or ESRD \n(dialysis) or SCr>3.3"
    } else if (outcome == "egfr40_2consec") {
      outcome_title = "ESRD or sustained (consecutive) \neGFR reduction of 40%"
    } else {
      outcome_title = dat_dict$long[which(dat_dict$short == outcome)]
    }
    p1 = ggsurvplot(km_fits_by_quart[[outcome]][[i]], data = df_sub,  censor = FALSE, 
                    title = paste0(outcome_title, ": \nQuartile ", i), 
                    conf.int = TRUE, risk.table = TRUE, fun = "event", 
                    legend.labs = c("Control", "Treatment"))
    
    # Create annotation for hazard ratio and 7-year RMST
    p1$plot = p1$plot + 
      annotate("label", x = 0, y = Inf, vjust = 1, hjust = 0, 
               size = 4, fontface = 2, 
               label = paste0("HR = ", 
                            est_with_ci_text(coxph_coeff_tabs[[outcome]][1+i,-1]), 
                            "\n7-year RMST difference\n   = ", 
                            est_with_ci_text(kfre_quart_rmst_boot_tabs[[outcome]][1+i,-1])))
    return(p1)
  })
  
  # Modify y-axis limits so that they are the same for all plots
  y_limits = sapply(all_plots, function(x) {
    layer_scales(x$plot)$y$range$range
  })
  for (i in 1:length(all_plots)) {
    all_plots[[i]]$plot = all_plots[[i]]$plot + 
      ylim(min(y_limits[1,]), max(y_limits[2,]))
  }
  
  # Collect all quartile plots for outcome
  png(paste0("figs_and_tabs/cum_inc_quart_", outcome, ".png"), 
       res = 100, width = 1200, height = 1200)
  arrange_ggsurvplots(all_plots, print = TRUE,
                      nrow = 2, ncol = 2)
  dev.off()
  arrange_ggsurvplots(all_plots, print = TRUE,
                      nrow = 2, ncol = 2)
}))

# Extract absolute risk for each arm, overall and by quartile
abs_risk_tabs = lapply(names(outcome_list), function(outcome){
  out = rbind(
    data.frame(Quartile = "Overall", 
               t(1 - summary(km_fits[[outcome]], times = horizon)$surv)
               ), 
    data.frame(Quartile = 1:4, 
               do.call(rbind, 
                       lapply(km_fits_by_quart[[outcome]], function(x){
                         t(1 - summary(x, times = horizon)$surv)
                         }))
               ))
  rownames(out) = NULL
  colnames(out) = c("Quartile", "Control", "Treatment")
  return(out)
})
names(abs_risk_tabs) = names(outcome_list)

# Function to calculate absolute risk reduction
calc_arr = function(km_fit, time, alpha = 0.05) {
  km_fit_summary = summary(km_fit, times = horizon)
  est = diff(km_fit_summary$surv)
  se = sqrt(sum(km_fit_summary$std.err^2))
  return(est + c(0, -1, 1) * qnorm(1 - alpha/2) * se)
}

# Calculate absolute risk reduction, overall and by quartile
arr_risk_tabs = lapply(names(outcome_list), function(outcome){
  out = rbind(
    data.frame(Quartile = "Overall", 
               t(calc_arr(km_fits[[outcome]], time = horizon))
               ), 
    data.frame(Quartile = 1:4, 
               do.call(rbind, 
                       lapply(km_fits_by_quart[[outcome]], function(x){
                         t(calc_arr(x, time = horizon))
                         }))
               ))
  rownames(out) = NULL
  colnames(out) = c("Quartile", "Est", "X2.5.", "X97.5.")
  return(out)
})
names(arr_risk_tabs) = names(outcome_list)

# Extract RMST for each arm, overall and by quartile
rmst_by_arm_tabs = lapply(names(outcome_list), function(outcome){
  out = rbind(
    data.frame(Quartile = "Overall", 
               Control = overall_rmst_fits[[outcome]]$RMST.arm0$rmst[[1]], 
               Treatment = overall_rmst_fits[[outcome]]$RMST.arm1$rmst[[1]]
               ), 
    data.frame(Quartile = 1:4, 
               do.call(rbind, 
                       lapply(kfre_quart_rmst_fits[[outcome]], function(x){
                         t(c(Control = x$RMST.arm0$rmst[[1]], 
                             Treatment = x$RMST.arm1$rmst[[1]]))
                         }))
               ))
  rownames(out) = NULL
  return(out)
})
names(rmst_by_arm_tabs) = names(outcome_list)

# Combine HR, absolute risks, ARR, RMSTs, and RMST difference into table
all_est_tab = lapply(names(outcome_list), function(outcome){
  out = data.frame(
    Quartile = c("Overall", 1:4), 
    t(sapply(1:5, function(i) {
      c(est_with_ci_text(coxph_coeff_tabs[[outcome]][i,-1]), 
        t(round(abs_risk_tabs[[outcome]][i,-1], 2)), 
        est_with_ci_text(arr_risk_tabs[[outcome]][i,-1]), 
        t(round(rmst_by_arm_tabs[[outcome]][i,-1], 2)), 
        est_with_ci_text(kfre_quart_rmst_boot_tabs[[outcome]][i,-1]))
    }))
  )
  rownames(out) = NULL
  colnames(out) = c("Quartile", "HR (95% CI)", 
                    "Event rate (control)", "Event rate (treatment)", 
                    "ARR (95% CI)", 
                    "RMST (control)", "RMST (treatment)", 
                    "RMST difference (95% CI)")
  return(out)
})
names(all_est_tab) = names(outcome_list)

# Format table for printing
est_tab = data.frame(
  Outcome = rep(names(all_est_tab), 
                each = nrow(all_est_tab[[1]])), 
  do.call(rbind, all_est_tab), 
  row.names = NULL
)
est_tab = merge(dat_dict[,c("short", "long")], est_tab, 
                  by.x = "short", by.y = "Outcome", all.y = TRUE, 
                  sort = FALSE)[,-1]
names(est_tab) = c("Outcome", "Quartile", "HR (95% CI)", 
                   "Event rate (control)", "Event rate (treatment)", 
                   "ARR (95% CI)", 
                   "RMST (control)", "RMST (treatment)", 
                   "RMST difference (95% CI))")
write.csv(est_tab, 
          file = "figs_and_tabs/est_tab.csv", row.names = FALSE)
est_tab
##                                                  Outcome Quartile
## 1                               Composite kidney outcome        1
## 2                               Composite kidney outcome        2
## 3                               Composite kidney outcome        3
## 4                               Composite kidney outcome  Overall
## 5                               Composite kidney outcome        4
## 6                       Development of macro-albuminuria  Overall
## 7                       Development of macro-albuminuria        1
## 8                       Development of macro-albuminuria        2
## 9                       Development of macro-albuminuria        4
## 10                      Development of macro-albuminuria        3
## 11           Renal failure or ESRD (dialysis) or SCr>3.3  Overall
## 12           Renal failure or ESRD (dialysis) or SCr>3.3        1
## 13           Renal failure or ESRD (dialysis) or SCr>3.3        2
## 14           Renal failure or ESRD (dialysis) or SCr>3.3        4
## 15           Renal failure or ESRD (dialysis) or SCr>3.3        3
## 16                                       All cause death  Overall
## 17                                       All cause death        1
## 18                                       All cause death        2
## 19                                       All cause death        3
## 20                                       All cause death        4
## 21                                  Cardiovascular death  Overall
## 22                                  Cardiovascular death        2
## 23                                  Cardiovascular death        3
## 24                                  Cardiovascular death        4
## 25                                  Cardiovascular death        1
## 26                       1st assisted hypoglycemic event  Overall
## 27                       1st assisted hypoglycemic event        2
## 28                       1st assisted hypoglycemic event        3
## 29                       1st assisted hypoglycemic event        4
## 30                       1st assisted hypoglycemic event        1
## 31 ESRD or sustained (consecutive) eGFR reduction of 40%  Overall
## 32 ESRD or sustained (consecutive) eGFR reduction of 40%        1
## 33 ESRD or sustained (consecutive) eGFR reduction of 40%        2
## 34 ESRD or sustained (consecutive) eGFR reduction of 40%        3
## 35 ESRD or sustained (consecutive) eGFR reduction of 40%        4
## 36                                                  <NA>  Overall
## 37                                                  <NA>        1
## 38                                                  <NA>        2
## 39                                                  <NA>        3
## 40                                                  <NA>        4
##          HR (95% CI) Event rate (control) Event rate (treatment)
## 1  0.93 (0.63, 1.37)                 0.09                   0.05
## 2  0.78 (0.59, 1.04)                 0.16                   0.12
## 3  0.73 (0.54, 0.97)                 0.15                    0.1
## 4  0.75 (0.65, 0.86)                 0.16                   0.11
## 5  0.66 (0.53, 0.83)                 0.26                   0.18
## 6  0.72 (0.61, 0.84)                 0.12                   0.09
## 7  0.73 (0.39, 1.39)                 0.05                   0.02
## 8  0.77 (0.56, 1.06)                 0.12                    0.1
## 9  0.62 (0.49, 0.79)                 0.23                   0.15
## 10 0.77 (0.56, 1.07)                 0.12                   0.08
## 11 0.93 (0.73, 1.18)                 0.05                   0.04
## 12 1.09 (0.68, 1.76)                 0.03                   0.03
## 13 0.89 (0.52, 1.51)                 0.04                   0.04
## 14 1.03 (0.69, 1.53)                 0.07                   0.06
## 15  0.65 (0.38, 1.1)                 0.04                   0.02
## 16   1.2 (1.04, 1.4)                 0.11                   0.12
## 17 1.37 (0.91, 2.07)                 0.04                   0.06
## 18 1.15 (0.84, 1.58)                 0.09                   0.08
## 19  0.95 (0.7, 1.29)                  0.1                    0.1
## 20  1.3 (1.03, 1.64)                 0.16                   0.18
## 21  1.3 (1.05, 1.62)                 0.05                   0.06
## 22 1.48 (0.89, 2.48)                 0.03                   0.04
## 23 0.85 (0.55, 1.31)                 0.06                   0.04
## 24  1.4 (1.01, 1.95)                 0.09                    0.1
## 25 1.75 (0.88, 3.47)                 0.01                   0.04
## 26  2.91 (2.56, 3.3)                 0.09                   0.22
## 27  3.9 (2.85, 5.35)                 0.08                   0.18
## 28 2.99 (2.32, 3.84)                 0.09                   0.24
## 29 2.41 (1.97, 2.94)                 0.16                    0.3
## 30 3.05 (2.22, 4.17)                 0.05                   0.17
## 31  1.02 (0.9, 1.15)                 0.25                    0.3
## 32  0.92 (0.7, 1.22)                 0.18                   0.47
## 33 1.03 (0.81, 1.31)                  0.3                   0.24
## 34 0.92 (0.73, 1.17)                 0.27                   0.21
## 35 1.16 (0.93, 1.46)                 0.19                   0.34
## 36 0.91 (0.81, 1.03)                 0.18                   0.16
## 37 1.08 (0.77, 1.51)                 0.11                    0.1
## 38 0.94 (0.73, 1.23)                 0.15                   0.13
## 39  0.8 (0.62, 1.01)                 0.18                   0.14
## 40 0.89 (0.73, 1.09)                 0.26                   0.23
##            ARR (95% CI) RMST (control) RMST (treatment)
## 1        0.04 (0, 0.09)        2465.89          2476.03
## 2     0.04 (-0.02, 0.1)        2375.33          2415.39
## 3         0.05 (0, 0.1)        2378.07          2428.15
## 4     0.05 (0.02, 0.07)        2365.05          2413.44
## 5     0.08 (0.03, 0.13)        2211.48          2326.26
## 6     0.04 (0.01, 0.06)        2408.31           2450.8
## 7    0.03 (-0.01, 0.08)        2519.65          2530.64
## 8    0.02 (-0.03, 0.08)        2419.04          2450.36
## 9     0.08 (0.03, 0.12)        2243.62          2360.11
## 10   0.03 (-0.01, 0.08)        2418.51          2451.97
## 11   0.01 (-0.01, 0.02)        2500.38          2505.81
## 12      0 (-0.02, 0.02)        2504.05          2499.95
## 13      0 (-0.04, 0.04)        2509.13          2516.66
## 14   0.01 (-0.02, 0.04)        2486.71           2487.7
## 15       0.02 (0, 0.04)        2503.06          2521.27
## 16  -0.01 (-0.03, 0.01)        2453.08          2429.44
## 17  -0.02 (-0.06, 0.01)        2499.84          2478.81
## 18   0.01 (-0.03, 0.04)        2460.21          2440.68
## 19      0 (-0.04, 0.04)        2443.66          2454.72
## 20  -0.02 (-0.07, 0.03)        2415.23          2358.53
## 21  -0.01 (-0.03, 0.01)        2505.78          2490.36
## 22  -0.01 (-0.02, 0.01)        2519.32          2502.86
## 23   0.01 (-0.02, 0.04)         2493.1          2504.98
## 24  -0.02 (-0.05, 0.02)         2478.8          2443.79
## 25  -0.03 (-0.06, 0.01)        2536.11          2518.29
## 26 -0.13 (-0.15, -0.11)        2427.33          2194.36
## 27  -0.1 (-0.15, -0.05)        2475.67          2259.61
## 28 -0.15 (-0.19, -0.11)        2419.58           2172.5
## 29 -0.15 (-0.19, -0.11)        2334.45          2046.06
## 30 -0.11 (-0.16, -0.07)        2477.46          2305.87
## 31    -0.05 (-0.2, 0.1)        2337.86             2336
## 32  -0.28 (-0.73, 0.16)        2383.63           2403.6
## 33   0.06 (-0.12, 0.23)        2360.03           2355.5
## 34   0.06 (-0.13, 0.24)         2348.4          2361.24
## 35   -0.16 (-0.4, 0.09)         2352.3          2322.99
## 36   0.02 (-0.01, 0.06)        2363.29          2376.76
## 37   0.01 (-0.07, 0.08)        2455.98          2446.84
## 38   0.02 (-0.04, 0.08)        2384.37          2390.91
## 39   0.04 (-0.01, 0.09)        2346.18          2387.56
## 40   0.03 (-0.02, 0.09)        2270.13           2293.4
##     RMST difference (95% CI))
## 1        10.14 (-23.9, 43.13)
## 2         40.07 (-4.65, 84.7)
## 3         50.08 (7.95, 92.69)
## 4          48.4 (25.34, 69.6)
## 5      114.77 (58.14, 176.41)
## 6        42.49 (22.52, 61.31)
## 7        10.99 (-8.12, 30.61)
## 8        31.32 (-4.86, 68.63)
## 9      116.49 (62.97, 173.51)
## 10       33.46 (-6.22, 71.98)
## 11        5.43 (-7.02, 16.95)
## 12       -4.1 (-30.36, 19.86)
## 13       7.53 (-14.43, 31.42)
## 14        0.99 (-27.6, 27.66)
## 15        18.22 (-3.14, 39.3)
## 16      -23.64 (-42.2, -6.64)
## 17       -21.03 (-47.7, 6.66)
## 18     -19.53 (-53.51, 12.56)
## 19      11.06 (-20.23, 42.14)
## 20     -56.7 (-100.19, -17.5)
## 21     -15.43 (-28.26, -3.57)
## 22      -16.46 (-40.91, 4.34)
## 23      11.88 (-12.68, 39.15)
## 24     -35.01 (-65.91, -5.41)
## 25       -17.82 (-37.08, 0.4)
## 26 -232.97 (-259.47, -208.03)
## 27 -216.06 (-263.96, -168.89)
## 28 -247.08 (-302.38, -192.26)
## 29 -288.38 (-346.09, -227.82)
## 30 -171.59 (-217.89, -127.53)
## 31      -1.86 (-25.85, 19.87)
## 32      19.97 (-25.65, 66.44)
## 33      -4.53 (-52.13, 41.24)
## 34      12.84 (-30.46, 56.44)
## 35     -29.32 (-75.96, 16.64)
## 36       13.46 (-10.1, 34.64)
## 37       -9.14 (-42.1, 26.45)
## 38       6.54 (-37.73, 48.36)
## 39        41.38 (0.72, 86.08)
## 40      23.27 (-33.49, 77.24)

Deciles

We repeat the HTE analysis using KFRE deciles to define the subgroups. The tables summarize the point estimates and 95% bootstrap CIs for 7-year RMST differences and normalized RMST differences by KFRE decile.

# kfre decile thresholds
kfre_dec_thresh = quantile(df$kfre5, seq(0.1, 0.9, by = 0.1))
# Create groups based on thresholds
df$kfre_dec = as.numeric(cut(df$kfre5, 
                             breaks = c(-Inf, as.numeric(kfre_dec_thresh), Inf),
                             labels = c(1:(length(kfre_dec_thresh)+1))))

# RMST difference by decile
kfre_dec_rmst = sapply(outcome_list, function(x){
  df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", 'kfre_dec')])
  if (x["status"] == "egfr40_2consec") {
    horizon = 2282 
  } else {
    horizon = 365*7
  }
  sapply(1:10, function(i){
    rmst2(df_sub[df_sub$kfre_dec==i,x["time"]], 
          df_sub[df_sub$kfre_dec==i,x["status"]], 
          df_sub[df_sub$kfre_dec==i,"glyarm"], 
          tau = horizon)$unadjusted.result[1,1]
  })
})
boot_kfre_dec_rmst = foreach::foreach(
  b = 1:B, 
  .packages = c("survRM2")
  ) %dopar% {
    # Set seed
    set.seed(b)
    # Create bootstrap sample
    boot_idx = sample(nrow(df), replace= TRUE)
    df_boot = df[boot_idx,]
    
    # Calculate overall RMST difference
    overall_rmst = sapply(outcome_list, function(x){
      df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "kfre_dec")])
      rmst2(df_sub[,x["time"]], 
            df_sub[,x["status"]], 
            df_sub[,"glyarm"], 
            tau = horizon)$unadjusted.result[1,1]
      })
    
    # Calculate RMST difference for each decile
    kfre_dec_rmst = sapply(outcome_list, function(x){
      df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "kfre_dec")])
      if (x["status"] == "egfr40_2consec") {
        horizon = 2282 
      } else {
        horizon = 365*7
      }
      sapply(1:10, function(i){
        rmst2(df_sub[df_sub$kfre_dec==i,x["time"]], 
              df_sub[df_sub$kfre_dec==i,x["status"]], 
              df_sub[df_sub$kfre_dec==i,"glyarm"], 
              tau = horizon)$unadjusted.result[1,1]
      })
    })
    return(list(overall_rmst = overall_rmst, 
                kfre_dec_rmst = kfre_dec_rmst))
  }

save(boot_kfre_dec_rmst, file = "boot_kfre_dec_rmst.rData")
# Not normalized
kfre_dec_rmst_boot_ci = lapply(1:10, function(i){
  dat = sapply(boot_kfre_dec_rmst, function(x){
    x$kfre_dec_rmst[i,]
  })
  apply(dat, 1, quantile, c(0.025, 0.975))
})

# Make table
kfre_dec_rmst_boot_tabs = lapply(names(outcome_list), function(outcome){
  out = rbind(
    data.frame(Decile = "Overall", 
               "Est" = overall_rmst[outcome], 
               t(kfre_overall_rmst_boot_ci[,outcome])), 
    data.frame(Decile = 1:10, 
               "Est" = kfre_dec_rmst[,outcome], 
               t(sapply(kfre_dec_rmst_boot_ci, function(x){x[,outcome]}))))
  rownames(out) = NULL
  return(out)
})
names(kfre_dec_rmst_boot_tabs) = names(outcome_list)

# Format table for printing
rmst_dec_tab = data.frame(
  Outcome = rep(names(kfre_dec_rmst_boot_tabs), 
                each = nrow(kfre_dec_rmst_boot_tabs[[1]])), 
  do.call(rbind, kfre_dec_rmst_boot_tabs), 
  row.names = NULL
)
rmst_dec_tab = merge(dat_dict[,c("short", "long")], rmst_dec_tab, 
                     by.x = "short", by.y = "Outcome", all.y = TRUE, 
                     sort = FALSE)[,-1]
names(rmst_dec_tab) = c("Outcome", "Decile", "Est", "2.5%", "97.5%")
rmst_dec_tab[,c("Est", "2.5%", "97.5%")] = 
  round(rmst_dec_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(rmst_dec_tab, 
          file = "figs_and_tabs/rmst_dec_tab.csv", row.names = FALSE)
rmst_dec_tab
##                                                  Outcome  Decile     Est
## 1                               Composite kidney outcome Overall   48.40
## 2                               Composite kidney outcome       3    8.71
## 3                               Composite kidney outcome       4   20.11
## 4                               Composite kidney outcome       1   35.50
## 5                               Composite kidney outcome       2  -23.52
## 6                               Composite kidney outcome       7   19.82
## 7                               Composite kidney outcome       8   76.14
## 8                               Composite kidney outcome       5   74.50
## 9                               Composite kidney outcome       6   40.66
## 10                              Composite kidney outcome       9  117.08
## 11                              Composite kidney outcome      10  166.55
## 12                      Development of macro-albuminuria       1   18.89
## 13                      Development of macro-albuminuria       2    0.97
## 14                      Development of macro-albuminuria Overall   42.49
## 15                      Development of macro-albuminuria       5   66.17
## 16                      Development of macro-albuminuria       6   21.69
## 17                      Development of macro-albuminuria       3    0.69
## 18                      Development of macro-albuminuria       4    8.81
## 19                      Development of macro-albuminuria       9  123.55
## 20                      Development of macro-albuminuria      10  166.73
## 21                      Development of macro-albuminuria       7   -2.21
## 22                      Development of macro-albuminuria       8   76.28
## 23           Renal failure or ESRD (dialysis) or SCr>3.3       3   12.64
## 24           Renal failure or ESRD (dialysis) or SCr>3.3       4   13.29
## 25           Renal failure or ESRD (dialysis) or SCr>3.3 Overall    5.43
## 26           Renal failure or ESRD (dialysis) or SCr>3.3       1   12.80
## 27           Renal failure or ESRD (dialysis) or SCr>3.3       2  -32.08
## 28           Renal failure or ESRD (dialysis) or SCr>3.3       7   13.57
## 29           Renal failure or ESRD (dialysis) or SCr>3.3       8   -5.17
## 30           Renal failure or ESRD (dialysis) or SCr>3.3       5    1.91
## 31           Renal failure or ESRD (dialysis) or SCr>3.3       6   24.34
## 32           Renal failure or ESRD (dialysis) or SCr>3.3       9    1.01
## 33           Renal failure or ESRD (dialysis) or SCr>3.3      10   10.81
## 34                                       All cause death       5  -38.36
## 35                                       All cause death       6   47.71
## 36                                       All cause death Overall  -23.64
## 37                                       All cause death       1   23.81
## 38                                       All cause death       2  -57.84
## 39                                       All cause death       3   -8.49
## 40                                       All cause death       4  -20.53
## 41                                       All cause death       9  -42.16
## 42                                       All cause death      10  -27.86
## 43                                       All cause death       7  -32.43
## 44                                       All cause death       8  -59.67
## 45                                  Cardiovascular death       7   -2.13
## 46                                  Cardiovascular death       8  -51.38
## 47                                  Cardiovascular death Overall  -15.43
## 48                                  Cardiovascular death       1    5.14
## 49                                  Cardiovascular death       2  -28.29
## 50                                  Cardiovascular death       3  -12.18
## 51                                  Cardiovascular death       4  -20.96
## 52                                  Cardiovascular death       5  -28.43
## 53                                  Cardiovascular death       6   31.45
## 54                                  Cardiovascular death       9  -23.39
## 55                                  Cardiovascular death      10  -15.24
## 56                       1st assisted hypoglycemic event       9 -255.22
## 57                       1st assisted hypoglycemic event      10 -322.17
## 58                       1st assisted hypoglycemic event Overall -232.97
## 59                       1st assisted hypoglycemic event       1 -156.51
## 60                       1st assisted hypoglycemic event       2 -179.27
## 61                       1st assisted hypoglycemic event       3 -146.58
## 62                       1st assisted hypoglycemic event       4 -246.15
## 63                       1st assisted hypoglycemic event       5 -244.96
## 64                       1st assisted hypoglycemic event       6 -222.18
## 65                       1st assisted hypoglycemic event       7 -223.98
## 66                       1st assisted hypoglycemic event       8 -322.58
## 67 ESRD or sustained (consecutive) eGFR reduction of 40% Overall   -1.86
## 68 ESRD or sustained (consecutive) eGFR reduction of 40%       1   20.17
## 69 ESRD or sustained (consecutive) eGFR reduction of 40%       2  -31.85
## 70 ESRD or sustained (consecutive) eGFR reduction of 40%       3   38.99
## 71 ESRD or sustained (consecutive) eGFR reduction of 40%       4   -1.02
## 72 ESRD or sustained (consecutive) eGFR reduction of 40%       5  -21.61
## 73 ESRD or sustained (consecutive) eGFR reduction of 40%       6   -2.60
## 74 ESRD or sustained (consecutive) eGFR reduction of 40%       7    7.88
## 75 ESRD or sustained (consecutive) eGFR reduction of 40%       8  -20.27
## 76 ESRD or sustained (consecutive) eGFR reduction of 40%       9   -6.69
## 77 ESRD or sustained (consecutive) eGFR reduction of 40%      10  -10.62
## 78                                                  <NA>       3    7.20
## 79                                                  <NA> Overall   13.46
## 80                                                  <NA>       1  -21.03
## 81                                                  <NA>       2   34.48
## 82                                                  <NA>       7    5.00
## 83                                                  <NA>       4  -11.74
## 84                                                  <NA>       5  -18.58
## 85                                                  <NA>       6   52.00
## 86                                                  <NA>       8   37.16
## 87                                                  <NA>       9  -14.76
## 88                                                  <NA>      10   78.95
##       2.5%   97.5%
## 1    25.34   69.60
## 2   -55.43   68.12
## 3   -59.20   94.31
## 4    -9.60   78.13
## 5   -71.09   24.33
## 6   -45.20   82.04
## 7    -0.95  150.41
## 8     4.43  139.17
## 9   -20.18  108.79
## 10   39.82  199.48
## 11   51.28  295.15
## 12   -2.94   43.91
## 13  -29.93   29.66
## 14   22.52   61.31
## 15    3.37  126.58
## 16  -36.22   82.96
## 17  -47.23   46.09
## 18  -53.19   68.30
## 19   48.33  201.58
## 20   56.66  287.16
## 21  -63.40   57.33
## 22   11.44  143.69
## 23  -30.94   53.54
## 24  -28.24   58.48
## 25   -7.02   16.95
## 26  -26.46   50.40
## 27  -67.33    3.42
## 28  -20.96   45.46
## 29  -46.12   37.26
## 30  -26.69   30.71
## 31   -8.94   60.53
## 32  -27.47   27.56
## 33  -40.83   64.28
## 34  -92.86    7.20
## 35   -2.11   93.33
## 36  -42.20   -6.64
## 37   -8.16   56.33
## 38 -112.47   -2.86
## 39  -56.91   39.78
## 40  -74.17   33.21
## 41  -97.85   11.41
## 42  -99.91   42.92
## 43  -84.23   17.94
## 44 -116.84    0.35
## 45  -37.81   38.23
## 46  -95.68   -3.82
## 47  -28.26   -3.57
## 48  -12.84   24.99
## 49  -67.66   12.53
## 50  -45.90   20.38
## 51  -58.87   16.21
## 52  -64.31    3.39
## 53   -4.53   67.83
## 54  -67.87   16.48
## 55  -69.36   39.69
## 56 -350.61 -162.96
## 57 -427.96 -213.76
## 58 -259.47 -208.03
## 59 -230.74  -94.45
## 60 -246.64 -110.25
## 61 -218.61  -74.88
## 62 -328.18 -166.06
## 63 -316.52 -170.72
## 64 -302.43 -141.34
## 65 -317.99 -138.66
## 66 -408.37 -237.44
## 67  -25.85   19.87
## 68  -25.69   66.34
## 69  -85.17   20.11
## 70  -23.29   93.70
## 71  -63.47   61.74
## 72  -73.14   30.99
## 73  -53.35   57.10
## 74  -51.01   68.09
## 75  -79.31   36.21
## 76  -60.38   51.59
## 77  -74.03   42.28
## 78  -58.05   67.32
## 79  -10.10   34.64
## 80  -66.53   25.80
## 81  -21.42   96.05
## 82  -69.94   83.43
## 83  -81.93   54.57
## 84  -94.39   47.63
## 85  -10.32  113.49
## 86  -41.26  116.84
## 87  -91.99   55.99
## 88  -10.99  171.19
# Normalized (subtract overall RMST difference)
kfre_dec_rmst_norm_boot_ci = lapply(1:10, function(i){
  dat = sapply(boot_kfre_dec_rmst, function(x){
    x$kfre_dec_rmst[i,] - x$overall_rmst
  })
  apply(dat, 1, quantile, c(0.025, 0.975))
})

# Make a table
kfre_dec_rmst_boot_norm_tabs = lapply(names(outcome_list), function(outcome){
  out = rbind(
    data.frame(Decile = 1:10, 
               "Est" = kfre_dec_rmst[,outcome] - overall_rmst[outcome], 
               t(sapply(kfre_dec_rmst_norm_boot_ci, function(x){x[,outcome]}))))
  rownames(out) = NULL
  return(out)
})
names(kfre_dec_rmst_boot_norm_tabs) = names(outcome_list)

# Format table for printing
rmst_norm_dec_tab = data.frame(
  Outcome = rep(names(kfre_dec_rmst_boot_norm_tabs), 
                each = nrow(kfre_dec_rmst_boot_norm_tabs[[1]])), 
  do.call(rbind, kfre_dec_rmst_boot_norm_tabs), 
  row.names = NULL
)
rmst_norm_dec_tab = merge(dat_dict[,c("short", "long")], rmst_norm_dec_tab, 
                          by.x = "short", by.y = "Outcome", all.y = TRUE, 
                          sort = FALSE)[,-1]
names(rmst_norm_dec_tab) = c("Outcome", "Decile", "Est", "2.5%", "97.5%")
rmst_norm_dec_tab[,c("Est", "2.5%", "97.5%")] = 
  round(rmst_norm_dec_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(rmst_norm_dec_tab, 
          file = "figs_and_tabs/rmst_norm_dec_tab.csv", row.names = FALSE)
rmst_norm_dec_tab
##                                                  Outcome Decile    Est    2.5%
## 1                               Composite kidney outcome      6  -7.74  -68.29
## 2                               Composite kidney outcome      7 -28.58  -94.28
## 3                               Composite kidney outcome      8  27.74  -41.41
## 4                               Composite kidney outcome      1 -12.90  -61.26
## 5                               Composite kidney outcome      2 -71.92 -118.99
## 6                               Composite kidney outcome      3 -39.69 -107.56
## 7                               Composite kidney outcome      4 -28.29 -102.44
## 8                               Composite kidney outcome      5  26.10  -41.20
## 9                               Composite kidney outcome     10 118.15   11.23
## 10                              Composite kidney outcome      9  68.68   -2.81
## 11                      Development of macro-albuminuria     10 124.24   18.35
## 12                      Development of macro-albuminuria      2 -41.52  -72.46
## 13                      Development of macro-albuminuria      9  81.06   12.61
## 14                      Development of macro-albuminuria      1 -23.60  -51.80
## 15                      Development of macro-albuminuria      6 -20.80  -77.63
## 16                      Development of macro-albuminuria      3 -41.80  -89.53
## 17                      Development of macro-albuminuria      4 -33.68  -90.49
## 18                      Development of macro-albuminuria      5  23.68  -35.90
## 19                      Development of macro-albuminuria      7 -44.70 -102.40
## 20                      Development of macro-albuminuria      8  33.80  -31.43
## 21           Renal failure or ESRD (dialysis) or SCr>3.3      1   7.36  -31.30
## 22           Renal failure or ESRD (dialysis) or SCr>3.3      3   7.21  -33.70
## 23           Renal failure or ESRD (dialysis) or SCr>3.3      4   7.86  -31.76
## 24           Renal failure or ESRD (dialysis) or SCr>3.3      5  -3.52  -31.38
## 25           Renal failure or ESRD (dialysis) or SCr>3.3      2 -37.51  -72.18
## 26           Renal failure or ESRD (dialysis) or SCr>3.3      7   8.14  -23.65
## 27           Renal failure or ESRD (dialysis) or SCr>3.3      8 -10.60  -48.28
## 28           Renal failure or ESRD (dialysis) or SCr>3.3      9  -4.42  -32.54
## 29           Renal failure or ESRD (dialysis) or SCr>3.3      6  18.91  -13.70
## 30           Renal failure or ESRD (dialysis) or SCr>3.3     10   5.38  -41.88
## 31                                       All cause death      2 -34.20  -83.89
## 32                                       All cause death      3  15.15  -33.06
## 33                                       All cause death      4   3.11  -47.17
## 34                                       All cause death      1  47.45   12.51
## 35                                       All cause death      6  71.35   25.75
## 36                                       All cause death      7  -8.79  -58.82
## 37                                       All cause death      8 -36.03  -91.25
## 38                                       All cause death      5 -14.72  -68.10
## 39                                       All cause death     10  -4.22  -69.41
## 40                                       All cause death      9 -18.52  -70.92
## 41                                  Cardiovascular death      6  46.88   11.00
## 42                                  Cardiovascular death      7  13.30  -23.09
## 43                                  Cardiovascular death      5 -13.00  -49.17
## 44                                  Cardiovascular death      1  20.57    0.95
## 45                                  Cardiovascular death      2 -12.87  -50.65
## 46                                  Cardiovascular death      3   3.25  -29.87
## 47                                  Cardiovascular death      4  -5.54  -41.36
## 48                                  Cardiovascular death      9  -7.96  -48.73
## 49                                  Cardiovascular death     10   0.19  -48.58
## 50                                  Cardiovascular death      8 -35.95  -76.10
## 51                       1st assisted hypoglycemic event     10 -89.20 -186.05
## 52                       1st assisted hypoglycemic event      8 -89.61 -169.00
## 53                       1st assisted hypoglycemic event      9 -22.24 -114.53
## 54                       1st assisted hypoglycemic event      1  76.46    4.79
## 55                       1st assisted hypoglycemic event      2  53.70   -9.64
## 56                       1st assisted hypoglycemic event      3  86.39   14.38
## 57                       1st assisted hypoglycemic event      4 -13.17  -87.19
## 58                       1st assisted hypoglycemic event      5 -11.98  -81.70
## 59                       1st assisted hypoglycemic event      6  10.79  -64.52
## 60                       1st assisted hypoglycemic event      7   9.00  -79.46
## 61 ESRD or sustained (consecutive) eGFR reduction of 40%      1  22.03  -25.22
## 62 ESRD or sustained (consecutive) eGFR reduction of 40%      2 -29.99  -78.54
## 63 ESRD or sustained (consecutive) eGFR reduction of 40%      3  40.85  -18.43
## 64 ESRD or sustained (consecutive) eGFR reduction of 40%      4   0.84  -59.93
## 65 ESRD or sustained (consecutive) eGFR reduction of 40%      5 -19.75  -70.09
## 66 ESRD or sustained (consecutive) eGFR reduction of 40%      6  -0.74  -49.54
## 67 ESRD or sustained (consecutive) eGFR reduction of 40%      7   9.74  -47.52
## 68 ESRD or sustained (consecutive) eGFR reduction of 40%      8 -18.41  -71.85
## 69 ESRD or sustained (consecutive) eGFR reduction of 40%      9  -4.83  -58.53
## 70 ESRD or sustained (consecutive) eGFR reduction of 40%     10  -8.76  -68.69
## 71                                                  <NA>      6  38.53  -22.03
## 72                                                  <NA>      7  -8.46  -82.39
## 73                                                  <NA>      8  23.69  -46.55
## 74                                                  <NA>      9 -28.23  -99.27
## 75                                                  <NA>     10  65.49  -18.51
## 76                                                  <NA>      2  21.02  -32.91
## 77                                                  <NA>      3  -6.26  -66.91
## 78                                                  <NA>      4 -25.21  -92.36
## 79                                                  <NA>      1 -34.49  -81.03
## 80                                                  <NA>      5 -32.04 -101.18
##     97.5%
## 1   56.03
## 2   32.81
## 3   98.84
## 4   32.11
## 5  -22.45
## 6   15.23
## 7   41.62
## 8   86.52
## 9  232.10
## 10 150.40
## 11 238.20
## 12  -7.63
## 13 155.53
## 14   2.42
## 15  36.38
## 16   4.41
## 17  25.05
## 18  79.58
## 19  12.97
## 20  97.58
## 21  43.90
## 22  46.82
## 23  50.35
## 24  26.20
## 25  -4.79
## 26  37.52
## 27  28.83
## 28  22.22
## 29  53.74
## 30  57.68
## 31  20.44
## 32  62.04
## 33  56.46
## 34  82.73
## 35 115.58
## 36  44.69
## 37  20.33
## 38  27.66
## 39  58.15
## 40  32.59
## 41  84.01
## 42  51.78
## 43  17.93
## 44  42.00
## 45  28.57
## 46  34.99
## 47  33.20
## 48  29.62
## 49  49.38
## 50  10.60
## 51  13.03
## 52  -8.63
## 53  64.96
## 54 139.23
## 55 121.70
## 56 157.80
## 57  63.94
## 58  61.59
## 59  87.80
## 60  99.37
## 61  66.89
## 62  21.12
## 63  96.75
## 64  56.84
## 65  34.63
## 66  57.24
## 67  69.21
## 68  37.15
## 69  49.27
## 70  50.02
## 71  95.67
## 72  66.51
## 73  97.65
## 74  38.38
## 75 149.57
## 76  80.61
## 77  51.53
## 78  37.03
## 79  10.34
## 80  29.27

RMST differences.

for (outcome in names(kfre_dec_rmst_boot_tabs)) {
  ate_tab = kfre_dec_rmst_boot_tabs[[outcome]][1,]
  p1 = kfre_dec_rmst_boot_tabs[[outcome]][-1,] %>% 
    mutate(x = 1:10) %>% 
    ggplot(aes(x = x, y = as.numeric(Est))) +
    geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) + 
    geom_hline(yintercept = ate_tab[,"Est"], color = "grey") + 
    annotate("rect", xmin = -Inf, xmax = Inf, 
             ymin = ate_tab[,"X2.5."], ymax = ate_tab[,"X97.5."], 
             alpha = 0.2) + 
    scale_x_continuous(breaks = 1:10, labels = 1:10) + 
    theme(legend.position = 'bottom', legend.box = 'vertical') + 
    xlab("Decile") + ylab("RMST difference") +
    ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
  print(p1)
  ggsave(paste0("figs_and_tabs/rmst_dec_", outcome, ".png"), 
         width = 6, height = 4)
}

Normalized RMST differences.

for (outcome in names(kfre_dec_rmst_boot_norm_tabs)) {
  p1 = kfre_dec_rmst_boot_norm_tabs[[outcome]] %>% 
    mutate(x = 1:10) %>% 
    ggplot(aes(x = x, y = as.numeric(Est))) +
    geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) + 
    scale_x_continuous(breaks = 1:10, labels = 1:10) + 
    theme(legend.position = 'bottom', legend.box = 'vertical') + 
    xlab("Decile") + ylab("Normalized RMST difference") + 
    geom_hline(yintercept = 0, color = "grey") + 
    ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
  print(p1)
  ggsave(paste0("figs_and_tabs/rmst_norm_dec_", outcome, ".png"), 
         width = 6, height = 4)
}

Jiang’s Diabetic Kidney Disease (DKD) score

We use the score from Jiang (2020) “Establishment and Validation of a Risk Prediction Model for Early Diabetic Kidney Disease Based on a Systematic Review and Meta-Analysis of 20 Cohorts” to define quartile subgroups. There are 1327 missing values, mostly due to missing smoking status.

# Diabetic retinopathy variable
baselinehistoryandphyisical = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/4-Data_Sets-CRFs/4a-CRF_Data_Sets/csv/f07_baselinehistoryphysicalexam.csv")
baselinehistoryandphyisical = baselinehistoryandphyisical %>% 
  mutate(retinopathy = case_when(retpathy==1 | lrtpathy==1 | rrtpathy==1 ~ 1, 
                                 TRUE ~ 0))

# Merge in diabetic retinopathy and years of diabetes
df_dkd =  merge(df, 
            baselinehistoryandphyisical %>% 
              select(MaskID, retinopathy), 
            by.x = "maskid", by.y = "MaskID", all.x = TRUE)

# Calculate dkd score from Jiang 2020
df_dkd = df_dkd %>% 
  mutate(dkd = case_when(baseline_age <= 49 ~ 0, 
                         baseline_age <= 59 ~ 3, 
                         baseline_age > 59 ~ 6) + 
           case_when(bmi < 25 ~ 0, 
                     bmi < 30 ~ 1.5, 
                     bmi >= 30  ~ 3) + 
           case_when(smokelif == FALSE ~ 0, 
                     smokelif == TRUE ~ 4) + 
           case_when(retinopathy == FALSE ~ 0, 
                     retinopathy == TRUE ~ 3) + 
           case_when(hba1c < 7 ~ 0, 
                     hba1c < 8 ~ 1.5, 
                     hba1c < 9 ~ 3, 
                     hba1c >= 9 ~ 4.5) + 
           case_when(sbp < 130 ~ 0, 
                     sbp < 140 ~ 2, 
                     sbp < 150 ~ 4, 
                     sbp >= 150 ~ 6) + 
           case_when(hdl < 1.3*18 ~ 0, 
                     hdl >= 1.3*18 ~ 2.5) + 
           case_when(trig < 1.7*18 ~ 0, 
                     trig >= 1.7*18 ~ 4)+ 
           case_when(uacr < 10 ~ 0, 
                     uacr < 20 ~ 2, 
                     uacr >= 20 ~ 4))

df_dkd %>% summarise_all(~ sum(is.na(.)))
##   maskid glyarm kfrs neph_composite neph_composite_fu Neph2 Neph2Days Neph3
## 1      0      0    0           1321              1321  1321      1321     0
##   Neph3Days ALLDEATH alldeath_fu CVDEATH cvdeath_fu hypoglycemia
## 1        17        0           0       0          0            0
##   hypoglycemia_fu cvd_primary cvd_primary_fu female baseline_age raceth sbp dbp
## 1               0           0              0      0            0      0  44  44
##   pulsepres chol trig vldl ldl hdl alt cpk fpg hba1c ualb ucreat uacr
## 1        44    2    2    2   2   1   0   1   0    14    0      0    0
##   ckd2021GFR screat bmi smokelif anti_htn chol_lowering diuretic insulin
## 1          0      0   6     1266        0             0        0       0
##   oral_DM raceth0 raceth1 raceth2 raceth3 days_from_baseline egfr40_2consec
## 1       0       0       0       0       0                  0              0
##   egfr40_2consec_fu kfre5 kfre_quarts kfre_dec retinopathy  dkd
## 1                17     0           0        0           0 1327
# Drop NAs
df_dkd = df_dkd %>% 
  drop_na(dkd)

# DKD quartile thresholds
dkd_quart_thresh = quantile(df_dkd$dkd, seq(0.25, 0.75, by = 0.25))
# Create groups based on thresholds
df_dkd$dkd_quarts = as.numeric(cut(df_dkd$dkd, 
                               breaks = c(-Inf, as.numeric(dkd_quart_thresh), Inf),
                               labels = c(1:(length(dkd_quart_thresh)+1))))

This score seems to be somewhat correlated with 5-year KFRE.

# Spearman's correlation
cor(df_dkd$kfre5, df_dkd$dkd, use = "complete.obs", method = "spearman")
## [1] 0.2941996
# Scatterplot
df_dkd %>% 
  ggplot(aes(x = kfre5, y = dkd)) + 
  geom_point() + 
  scale_x_log10() + 
  xlab("5-year predicted risk by KFRE (log)") + ylab("Jiang's DKD score")

ggsave("figs_and_tabs/dkd_vs_kfre.png", width = 6, height = 5)

Number of observations overall and in each arm, restricting to those with non-missing covariates needed to calculate the DKD score.

c(table(df_dkd$glyarm), Overall = nrow(df_dkd))
##       0       1 Overall 
##    4231    4219    8450

Number of non-missing observations for each outcome.

t(sapply(outcome_list, function(x) {
  my_df = na.omit(df_dkd[,c(x["status"], x["time"], "glyarm")])
  c(table(my_df$glyarm), Overall = nrow(my_df))
}))
##                   0    1 Overall
## neph_composite 3708 3649    7357
## Neph2          3708 3649    7357
## Neph3          4224 4210    8434
## ALLDEATH       4231 4219    8450
## CVDEATH        4231 4219    8450
## hypoglycemia   4231 4219    8450
## cvd_primary    4231 4219    8450
## egfr40_2consec 4224 4210    8434

We estimate the 7-year RMST differences for each subgroup defined by the DKD quartiles.

# Overall
dkd_overall_rmst = sapply(outcome_list, function(x){
  df_sub = na.omit(df_dkd[,c(x["time"], x["status"], "glyarm", "dkd_quarts")])
  rmst2(df_sub[,x["time"]], 
        df_sub[,x["status"]], 
        df_sub[,"glyarm"], 
        tau = horizon)$unadjusted.result[1,1]
})
# By quartile
dkd_quart_rmst = sapply(outcome_list, function(x){
  df_sub = na.omit(df_dkd[,c(x["time"], x["status"], "glyarm", "dkd_quarts")])
  if (x["status"] == "egfr40_2consec") {
    horizon = 2548 
  } else {
    horizon = 365*7
  }
  sapply(1:4, function(i){
    rmst2(df_sub[df_sub$dkd_quarts==i,x["time"]], 
          df_sub[df_sub$dkd_quarts==i,x["status"]], 
          df_sub[df_sub$dkd_quarts==i,"glyarm"], 
          tau = horizon)$unadjusted.result[1,1]
  })
})

We use 1000 bootstrap samples to obtain confidence intervals for the RMST differences. Tables of the RMST differences and normalized RMST differences (where the overall RMST difference is subtracted from each quartile’s estime) with 95% bootstrap CIs are shown below.

boot_dkd_rmst = foreach::foreach(
  b = 1:B, 
  .packages = c("survRM2")
  ) %dopar% {
    # Set seed
    set.seed(b)
    # Create bootstrap sample
    boot_idx = sample(nrow(df), replace= TRUE)
    df_boot = df_dkd[boot_idx,]
    
    # Calculate overall RMST difference
    overall_rmst = sapply(outcome_list, function(x){
      df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "dkd_quarts")])
      rmst2(df_sub[,x["time"]], 
            df_sub[,x["status"]], 
            df_sub[,"glyarm"], 
            tau = horizon)$unadjusted.result[1,1]
      })
    
    # Calculate RMST difference for each quartile
    dkd_quart_rmst = sapply(outcome_list, function(x){
      df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "dkd_quarts")])
      if (x["status"] == "egfr40_2consec") {
        horizon = 2548 
      } else {
        horizon = 365*7
      }
      sapply(1:4, function(i){
        rmst2(df_sub[df_sub$dkd_quarts==i,x["time"]], 
              df_sub[df_sub$dkd_quarts==i,x["status"]], 
              df_sub[df_sub$dkd_quarts==i,"glyarm"], 
              tau = horizon)$unadjusted.result[1,1]
      })
    })
    return(list(overall_rmst = overall_rmst, 
                dkd_quart_rmst = dkd_quart_rmst))
  }

save(boot_dkd_rmst, file = "boot_dkd_rmst.rData")
# Overall
dkd_overall_rmst_boot_ci = apply(
  sapply(boot_dkd_rmst, function(x){
    x$overall_rmst
  }), 1, quantile, c(0.025, 0.975))

# Not normalized
dkd_quart_rmst_boot_ci = lapply(1:4, function(i){
  dat = sapply(boot_dkd_rmst, function(x){
    x$dkd_quart_rmst[i,]
  })
  apply(dat, 1, quantile, c(0.025, 0.975))
})

# Make tables
dkd_quart_rmst_boot_tabs = lapply(names(outcome_list), function(outcome){
  out = rbind(
    data.frame(Quartile = "Overall", 
               "Est" = dkd_overall_rmst[outcome], 
               t(dkd_overall_rmst_boot_ci[,outcome])), 
    data.frame(Quartile = 1:4, 
               "Est" = dkd_quart_rmst[,outcome], 
               t(sapply(dkd_quart_rmst_boot_ci, function(x){x[,outcome]}))))
  rownames(out) = NULL
  return(out)
})
names(dkd_quart_rmst_boot_tabs) = names(outcome_list)

# Format table for printing
dkd_rmst_tab = data.frame(
  Outcome = rep(names(dkd_quart_rmst_boot_tabs), 
                each = nrow(dkd_quart_rmst_boot_tabs[[1]])), 
  do.call(rbind, dkd_quart_rmst_boot_tabs), 
  row.names = NULL
)
dkd_rmst_tab = merge(dat_dict[,c("short", "long")], dkd_rmst_tab, 
                     by.x = "short", by.y = "Outcome", all.y = TRUE, 
                     sort = FALSE)[,-1]
names(dkd_rmst_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
dkd_rmst_tab[,c("Est", "2.5%", "97.5%")] = 
  round(dkd_rmst_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(dkd_rmst_tab, 
          file = "figs_and_tabs/dkd_rmst_tab.csv", row.names = FALSE)
dkd_rmst_tab
##                                                  Outcome Quartile     Est
## 1                               Composite kidney outcome        1   -1.45
## 2                               Composite kidney outcome        2   85.18
## 3                               Composite kidney outcome        3   42.25
## 4                               Composite kidney outcome  Overall   44.87
## 5                               Composite kidney outcome        4   66.30
## 6                       Development of macro-albuminuria  Overall   39.29
## 7                       Development of macro-albuminuria        1   13.15
## 8                       Development of macro-albuminuria        2   59.52
## 9                       Development of macro-albuminuria        4   74.78
## 10                      Development of macro-albuminuria        3   26.68
## 11           Renal failure or ESRD (dialysis) or SCr>3.3  Overall    2.02
## 12           Renal failure or ESRD (dialysis) or SCr>3.3        1  -18.98
## 13           Renal failure or ESRD (dialysis) or SCr>3.3        2   16.80
## 14           Renal failure or ESRD (dialysis) or SCr>3.3        4  -14.72
## 15           Renal failure or ESRD (dialysis) or SCr>3.3        3   23.12
## 16                                       All cause death  Overall  -24.05
## 17                                       All cause death        1  -21.62
## 18                                       All cause death        2  -13.31
## 19                                       All cause death        3  -30.15
## 20                                       All cause death        4  -31.30
## 21                                  Cardiovascular death  Overall  -14.34
## 22                                  Cardiovascular death        2  -13.14
## 23                                  Cardiovascular death        3   -7.32
## 24                                  Cardiovascular death        4  -19.76
## 25                                  Cardiovascular death        1  -17.02
## 26                       1st assisted hypoglycemic event  Overall -243.40
## 27                       1st assisted hypoglycemic event        2 -215.30
## 28                       1st assisted hypoglycemic event        3 -265.54
## 29                       1st assisted hypoglycemic event        4 -318.35
## 30                       1st assisted hypoglycemic event        1 -189.00
## 31 ESRD or sustained (consecutive) eGFR reduction of 40%  Overall   -6.39
## 32 ESRD or sustained (consecutive) eGFR reduction of 40%        1  -22.44
## 33 ESRD or sustained (consecutive) eGFR reduction of 40%        2   18.58
## 34 ESRD or sustained (consecutive) eGFR reduction of 40%        3   35.63
## 35 ESRD or sustained (consecutive) eGFR reduction of 40%        4  -61.36
## 36                                                  <NA>  Overall   14.79
## 37                                                  <NA>        1    2.73
## 38                                                  <NA>        2   21.86
## 39                                                  <NA>        3   -1.66
## 40                                                  <NA>        4   44.36
##       2.5%   97.5%
## 1   -31.25   31.69
## 2    46.13  126.42
## 3    -6.62   98.21
## 4    22.40   70.88
## 5     3.02  135.12
## 6    20.68   61.41
## 7    -9.41   36.59
## 8    26.35   92.45
## 9    14.07  137.76
## 10  -15.88   76.56
## 11  -10.29   15.15
## 12  -41.51    3.25
## 13   -9.23   42.69
## 14  -43.88   16.34
## 15   -3.98   51.16
## 16  -42.16   -5.43
## 17  -46.72    4.35
## 18  -46.95   19.37
## 19  -69.00    7.03
## 20  -78.55   17.71
## 21  -27.16   -0.42
## 22  -41.48   10.65
## 23  -34.33   19.71
## 24  -55.01   17.40
## 25  -35.28    1.49
## 26 -273.12 -215.90
## 27 -264.70 -167.93
## 28 -320.57 -206.13
## 29 -386.64 -249.73
## 30 -237.16 -144.80
## 31  -29.37   18.53
## 32  -62.25   18.13
## 33  -22.70   62.02
## 34  -14.47   87.44
## 35 -122.34   -3.70
## 36   -8.20   38.68
## 37  -33.68   40.86
## 38  -21.75   65.14
## 39  -52.46   46.67
## 40  -18.26  102.03
# Normalized (subtract overall RMST difference)
dkd_quart_rmst_norm_boot_ci = lapply(1:4, function(i){
  dat = sapply(boot_dkd_rmst, function(x){
    x$dkd_quart_rmst[i,] - x$overall_rmst
  })
  apply(dat, 1, quantile, c(0.025, 0.975))
})

# Make tables
dkd_quart_rmst_boot_norm_tabs = lapply(names(outcome_list), function(outcome){
  data.frame(Quartile = 1:4, 
             "Est" = dkd_quart_rmst[,outcome] - dkd_overall_rmst[outcome], 
             t(sapply(dkd_quart_rmst_norm_boot_ci, function(x){x[,outcome]})))
})
names(dkd_quart_rmst_boot_norm_tabs) = names(outcome_list)

# Format table for printing
dkd_rmst_norm_tab = data.frame(
  Outcome = rep(names(dkd_quart_rmst_boot_norm_tabs), 
                each = nrow(dkd_quart_rmst_boot_norm_tabs[[1]])), 
  do.call(rbind, dkd_quart_rmst_boot_norm_tabs), 
  row.names = NULL
)
dkd_rmst_norm_tab = merge(dat_dict[,c("short", "long")], dkd_rmst_norm_tab, 
                          by.x = "short", by.y = "Outcome", all.y = TRUE, 
                          sort = FALSE)[,-1]
names(dkd_rmst_norm_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
dkd_rmst_norm_tab[,c("Est", "2.5%", "97.5%")] = 
  round(dkd_rmst_norm_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(dkd_rmst_norm_tab, 
          file = "figs_and_tabs/dkd_rmst_norm_tab.csv", row.names = FALSE)
dkd_rmst_norm_tab
##                                                  Outcome Quartile    Est
## 1                               Composite kidney outcome        1 -46.32
## 2                               Composite kidney outcome        2  40.30
## 3                               Composite kidney outcome        3  -2.62
## 4                               Composite kidney outcome        4  21.42
## 5                       Development of macro-albuminuria        1 -26.14
## 6                       Development of macro-albuminuria        2  20.23
## 7                       Development of macro-albuminuria        3 -12.61
## 8                       Development of macro-albuminuria        4  35.49
## 9            Renal failure or ESRD (dialysis) or SCr>3.3        1 -21.00
## 10           Renal failure or ESRD (dialysis) or SCr>3.3        2  14.78
## 11           Renal failure or ESRD (dialysis) or SCr>3.3        3  21.11
## 12           Renal failure or ESRD (dialysis) or SCr>3.3        4 -16.74
## 13                                       All cause death        1   2.43
## 14                                       All cause death        2  10.74
## 15                                       All cause death        3  -6.11
## 16                                       All cause death        4  -7.25
## 17                                  Cardiovascular death        1  -2.68
## 18                                  Cardiovascular death        2   1.20
## 19                                  Cardiovascular death        3   7.02
## 20                                  Cardiovascular death        4  -5.43
## 21                       1st assisted hypoglycemic event        1  54.40
## 22                       1st assisted hypoglycemic event        2  28.11
## 23                       1st assisted hypoglycemic event        3 -22.14
## 24                       1st assisted hypoglycemic event        4 -74.95
## 25 ESRD or sustained (consecutive) eGFR reduction of 40%        1 -16.05
## 26 ESRD or sustained (consecutive) eGFR reduction of 40%        2  24.97
## 27 ESRD or sustained (consecutive) eGFR reduction of 40%        3  42.02
## 28 ESRD or sustained (consecutive) eGFR reduction of 40%        4 -54.97
## 29                                                  <NA>        1 -12.06
## 30                                                  <NA>        2   7.07
## 31                                                  <NA>        3 -16.45
## 32                                                  <NA>        4  29.57
##       2.5%  97.5%
## 1   -79.03 -17.43
## 2     3.42  78.62
## 3   -43.95  42.97
## 4   -33.15  77.88
## 5   -53.10  -1.76
## 6    -9.15  50.15
## 7   -49.53  27.74
## 8   -15.36  88.31
## 9   -41.77  -0.67
## 10   -8.62  36.78
## 11   -1.54  44.95
## 12  -41.99   9.48
## 13  -24.89  28.70
## 14  -18.88  37.15
## 15  -38.22  25.78
## 16  -46.67  34.86
## 17  -21.48  17.29
## 18  -22.85  21.44
## 19  -18.10  29.57
## 20  -35.55  27.46
## 21    9.45  94.72
## 22  -14.99  70.82
## 23  -70.53  25.73
## 24 -131.28 -16.13
## 25  -53.80  19.94
## 26  -13.72  59.54
## 27   -1.53  86.23
## 28 -108.37  -5.96
## 29  -46.23  22.40
## 30  -29.95  42.04
## 31  -59.41  23.36
## 32  -21.98  78.94

These are plots of the RMST differences by outcome, with a horizontal reference line drawn at the overall/ATE point estimate. Shading denotes the ATE CI.

for (outcome in names(dkd_quart_rmst_boot_tabs)) {
  ate_tab = dkd_quart_rmst_boot_tabs[[outcome]][1,]
  p1 = dkd_quart_rmst_boot_tabs[[outcome]][-1,] %>% 
    mutate(x = 1:4) %>% 
    ggplot(aes(x = x, y = as.numeric(Est))) +
    geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) + 
    geom_hline(yintercept = ate_tab[,"Est"], color = "grey") + 
    annotate("rect", xmin = -Inf, xmax = Inf, 
             ymin = ate_tab[,"X2.5."], ymax = ate_tab[,"X97.5."], 
             alpha = 0.2) + 
    scale_x_continuous(breaks = 1:4, labels = 1:4) + 
    theme(legend.position = 'bottom', legend.box = 'vertical') + 
    xlab("Quartile") + ylab("RMST difference") +
    ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
  print(p1)
  ggsave(paste0("figs_and_tabs/dkd_rmst_", outcome, ".png"), 
         width = 6, height = 4)
}

Normalized RMST differences.

for (outcome in names(dkd_quart_rmst_boot_norm_tabs)) {
  p1 = dkd_quart_rmst_boot_norm_tabs[[outcome]] %>% 
    mutate(x = 1:4) %>% 
    ggplot(aes(x = x, y = as.numeric(Est))) +
    geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) + 
    scale_x_continuous(breaks = 1:4, labels = 1:4) + 
    theme(legend.position = 'bottom', legend.box = 'vertical') + 
    xlab("Quartile") + ylab("Normalized RMST difference") + 
    geom_hline(yintercept = 0, color = "grey") + 
    ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
  print(p1)
  ggsave(paste0("figs_and_tabs/dkd_rmst_norm_", outcome, ".png"), 
         width = 6, height = 4)
}